T6-A181  7M 


UNCLASSIFIED 


FURTHER  DEVELOPMENT  OF  A  ONE -DIMENSIONAL  UNSTEAOV  EULER 
CODE  FOR  HAVE  ROTOR  APPLICATIONS^)  NAVAL  POSTGRADUATE 
SCHOOL  NONTEREV  CA  D  T  JOHNSTON  MAR  87  NP567-87-88 2 


MICROCOPY  RESOLUTION  TEST  CHART 


AD-A181  750 


#  *  ^ 


NPS67-87-002 


NAVAL  POSTGRADUATE  SCHOOL 


Monterey,  California 


•  r  '  •  r  C'f- r» 


OTIC 

1L.ECTE 
m  2  9  ^7 


JD 


THESIS 


FURTHER  DEVELOPMENT  OF  A 
ONE-DIMENSIONAL  UNSTEADY  EULER  CODE 
FOR  WAVE  ROTOR  APPLICATIONS 


by 


David  T.  Johnston 


March  1987 


Thesis  Advisor: 


Ray  P .  Shreeve 


Approved  for  public  release;  distribution  is  unlimited 

Prepared  for:  Naval  Air  Systems  Command 

Washington,  DC  20361 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


Rear  Admiral  Robert  C.  Austin  David  A.  Schrady 

Superintendent  Provost 


This  thesis  was  prepared  in  conjunction  with  research 
sponsored  in  part  by  the  Naval  Air  Systems  Command, 
Washington,  DC,  under  Air  Task  #A931931E/186A/7R024-03-001. 

Reproduction  of  all  or  part  of  this  work  is  authorized. 


i  Released  by: 


GORDON  E.  SCHACHER 
Dean  of  Science  and 


Engineering 


UNCLASSIFIED 


u  report  security  classification 

UNCLASSIFIED 


it  SECURlT t  C. ASSlFlCATlON  AUTHORITY 


//  7. 


5 


REPORT  DOCUMENTATION  PAGE 


lb  RESTRICTIVE  MARKINGS 


Zb  OEClASSiFiCATiON/ DOWNGRADING  SCHEDULE 


4  PERFORMING  ORGANIZATION  REPORT  NUMSER(S) 

NPS67-87-002 


64  NAME  OF  PERFORMING  ORGANIZATION 


)  DISTRIBUTION/  AVAILABILITY  OF  REPORT 

Approved  for  public  release; 
distribution  is  unlimited 


S  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

NPS67-87-002 


6b  OFFICE  SYMBOL  >4  NAME  OF  MONITORING  ORGANIZATION 

(It  400/«4b/*) 


Naval  Postgraduate  School  Code  67 


6 <  AOORESS  (Cry.  Sfvr#.  trtd  HP Cod*) 


Naval  Postgraduate  School 


?b  ADDRESS  (Ciry.  SUf*.  and  Z/P  Cod*) 


Monterey,  California  93943-5000 


Monterey,  California  93943-5000 


8b  OFFICE  SYMBOL  9  PROCUREMENT  INSTRUMENT  lOEN Tif iCA TION  NUMBER 
(If  sppliCible) 

AIR  931E  NO 0 0 1 9 8 7V7R7 R04 8 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 
ELEMENT  NO 


CT 

task 

WORK  jNi r 

V 

NO 

ACCESSION  NO 

24 

03 

001 

8j  NAME  OF  FUNOING/ SPONSORING  1 8b  OFFICE  SYMI 

ORGANIZATION  I  (It  tpplidble 

Naval  Air  Systems  Command!  AIR  9  31E 


8c  AOORESS  (Ciry.  Sf4f*.  end  HPC od*J 


61153N 


M  TiTlE  (Include  Security  Cltaihcetion) 

FURTHER  DEVELOPMENT  OF  A  ONE- DIMENSIONAL  UNSTEADY  EULER  CODE  FOR  WAVE 
ROTOR  APPLICATIONS 


tZ  PERSONAL  auThOR(S) 

Johnston,  David  T.  / 


'34  type  of  REPORT 

Master's  Thesis 


■  6  supplementary  notation 


14  DATE  OF  REPORT  (Year  Month  Oey)  IS  PAGE  COuNT 

1987,  March  204 


•  Y 

COSATI  COOES  i 

F  ElO 

GROUP 

SUBGROUP  | 

mmmmm 

18  VbVJECL  JipMS  1 Continue  on  reverie  if  necetssry jpd  identify  Qy  Slock  number) 

■QAZ4 (v;  Wave  Rotor;  Uasiefray  Euler  Code; 
Computational  Fluid  Dynamics 


m 


PACT  Continue  on  reverie  <t  n*<*U4«y  jno  identity  by  bloc*  numtoerl 

The  EULERl  Fortran  program  for  computing  one-dimensional  unsteady  flow 
based  on  the  QAZlD  method  of  Verhoff  was  extended  to  provide  tracking  and 
correcting  of  discontinuities,  flow  to  right  or  left  and  open  and  closed 
end  boundary  conditions.  The  program  was  run  on  four  shock  tube  test 
cases  to  verify  accuracy  and  range  of  capability  of  the  revised  E1DV2 
code.  Further  extensions  and  test  verifications  are  recommended  so  that 
the  code  is  suitable  for  wave  rotor  applications,  ; 


10  0  S'RiBUTlON/ AVAILABILITY  of  ABSTRACT 

rjtuNCLASSlFlEOAjNLiMlTED  □  SAME  AS  RPT  □  OTIC  USERS 

Z1  ABSTRACT  security  CLASSIFICATION 

Unclassified 

<V4  NAME  OF  RESPONSIBLE  iNO'ViOUAL 

Prof.  Raymond  P.  Shreeve  . 

ZZb  TELE Ph 

(408) 

iONE  (It 

_ §4 

iclude  Ares  Code)  j 

6-2593 _ 

mmmmm 

00  FORM  1473,84 mar 


83  APR  bdition  m»y  b*  utcd  until  *ih»uit«d 
All  other  cd'hont  4tV  obiolftt 
1 


SECURITY  CLASSIFICATION  of  Tm$  page 

UNCLASSIFIED 


i 


Approved  for  public  release;  distribution  is  unlimited 


Further  Development  of  a 
One-Dimensional  Unsteady  Euler  Code 
for  Wave  Rotor  Applications 

by 

David  T.  Johnston 
Lieutenant,  United  States  Navy 
B.S.,  Cornell  University,  1979 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  AERONAUTICAL  ENGINEERING 


from  the 

NAVAL  POSTGRADUATE  SCHOOL 
March  1987 


JZU/LrJL^ 

G.E.  Schacher, 

Dean  of  Science  and  Engineering 


unsteady  flow  based  on  the  QAZ1D  method  of  Verhoff  was 
extended  to  provide  tracking  and  correcting  of  discontinui¬ 
ties,  flow  to  right  or  left  and  open  and  closed  end  boundary 


conditions.  The  program  was  run  on  four  shock  tube  test 
cases  to  verify  accuracy  and  range  of  capability  of  the 
revised  E1DV2  code.  Further  extensions  and  test  verifica¬ 
tions  are  recommended  so  that  the  code  is  suitable  for  wave 
rotor  applications. 

Acce:>ion  for 

MIS  CffA&l 
DTIC  TAB 

U  .anno.i'vvd 
J’.iatdic  jI.um 


f~ 


I.  INTRODUCTION -  18 

II.  ANALYTICAL  AND  COMPUTATIONAL  APPROACH  -  20 

A.  QAZ1D  DESCRIPTION  -  20 

B.  CHARACTERISTIC  SOLUTION  -  23 

C.  DISCONTINUITIES -  28 

D.  BOUNDARY  CONDITIONS  -  44 

III.  FORTRAN  PROGRAM  E1DV2  -  51 

A.  GENERAL  DESCRIPTION  -  51 

B.  THE  MAIN  PROGRAM -  52 

C.  THE  SUBROUTINE  PROGRAMS -  56 

IV.  RESULTS -  76 

A.  TEST  CASE  1 -  77 

B.  TEST  CASE  2  -  83 

C.  TEST  CASE  3  -  36 

D.  TEST  CASE  4  -  36 

V.  DISCUSSION -  92 

A.  RESULTS  OF  TEST  CASES -  92 

B.  CURRENT  LIMITATIONS  OF  THE  E1DV2  CODE  -  96 

VI.  CONCLUSIONS -  99 

APPENDIX  A:  DERIVATION  OF  EQUATIONS  -  101 

APPENDIX  B:  OPERATION  OF  ED1V2  ON  THE  NPS  VM/CMS 

SYSTEM -  115 

APPENDIX  C:  FLOWCHARTS  -  122 

APPENDIX  D:  E1DV2  FORTRAN  LISTING  -  143 


LIST  OF  REFERENCES -  200 

INITIAL  DISTRIBUTION  LIST -  201 


5 


SUMMARY  OF  THE  EULER  EQUATIONS  - 

SHOCK  EQUATIONS  FOR  HIGH  PRESSURE  ON  LEFT 
SHOCK  EQUATIONS  FOR  HIGH  PRESSURE  ON  RIGHT 

SUBROUTINES  - 

VALUES  OF  PARAMETERS  SHOCK  AND  CNTACT  - 

TERMINAL  AND  OUTPUT  - 


EDITABLE  VARIABLES 


2 . 1  The  Natural  Coordinate  System  - 

2.2  Characteristic  Curve  within  Computational  Mesh  — 

2.3  Computational  Grid  for  the  QAZ1D  Method  - 


2.4  Shock  Wave  with  High  Pressure  on  Left  - 

2.5  Q  Extended  Riemann  Variable  Change  with 

Mach  Number  - 

2 . 6  Shock  Wave  with  High  Pressure  on  Right  - 

2 . 7  Extended  Riemann  Variable  Change  with  Mach 

Number,  High  Pressure  on  Right  - 

2.8  Behavior  of  Physical  Properties  through  a 

Contact  Surface  - 

2.9  Modified  Entropy  Distribution  for  Shock  and 

Contact  Surface  Traveling  Right  - 

2.10  Modified  Entropy  Distribution  for  Shock  and 

Contact  Surface  Traveling  Left  - 

2.11  Left  Boundary  Computational  Grid  - 

2.12  Right  Boundary  Computational  Grid  - 

2.13  Shock  Reflecting  at  Solid  Boundary  - 

3.1  Overall  Algorithm  for  QAZ1D  Method  - 

3.2  Interval  and  Node  Description  used  in  "SWEEP"  - 

3.3  Schematic  of  SHOCK  =  331  and  CNTACT  =  232  - 

3.4  Schematic  of  SHOCK  =  100  and  CNTACT  =  321  - 

3.5  Schematic  of  SHOCK  =  221  and  CNTACT  =  222  - 

3.6  Schematic  of  SHOCK  =  231  and  CNTACT  =  322  - 

3.7  Computational  Method  for  "CONDI"  Routine  - 

3.8  Computational  Method  for  "C0ND2"  Routine  - 


3.9 


Computational  Method  for  "C0ND3"  Routine 


67 


3.10  Example  Condition  for  "COND5"  Routine - 

3.11  Example  of  "C0ND4"  Routine  Situation  - 

3.12  Example  of  "CONDe”  Situation  - 

3 . 13  Discontinuity  Location  within  an  Interval 

4 . 1  Shock  Tube  at  t  =  0  and  t  =  t  - 

4.2  Test  Case  1  (J  =  1  to  J  =  55)  - 

4.3  Test  Case  1  (J  =  56  to  J  =  109)  - 

4.4  Exact  and  Computed  Density  Distributions 

(Test  Case  1)  - 

4.5  Location  of  Discontinuities  versus  Time 

(Test  Case  1)  - 

4.6  Test  Case  1,  High  Pressure  on  Right  - 

4 . 7  Test  Case  2 ,  High  Pressure  on  Left  - 

4 . 8  Test  Case  2 ,  High  Pressure  on  Right  - 

4.9  Test  Case  3,  High  Pressure  on  Left  - 

4 . 10  Test  Case  3 ,  High  Pressure  on  Right  - 

4.11  Test  Case  4,  High  Pressure  on  Left  - 

4.12  Test  Case  4,  High  Pressure  on  Right  - 

A. 1  Shock  Wave  with  High  Pressure  on  Right  - 

A. 2  Contact  Surface  Traveling  Right,  Subscript 
Notation  - 

A. 3  Contact  Surface  Traveling  Left,  Subscript 
Notation  - 

C.l  Main  Program  Flowchart  - 

C.2  "DBURST"  Subroutine  Flowchart  - 

C.3  "TRAK"  Subroutine  Flowchart  - 


68 

69 

70 
73 
76 

78 

79 

80 

81 

82 

84 

85 
37 
88 

90 

91 
102 


112 
123 
12  5 
127 


C.  4 


"SWEEP"  Subroutine  Flowchart 


C . 5  General  "CONDI, 2, 3  and  5"  Subroutine  Flowchart  —  134 

C. 6  "C0ND4"  Subroutine  Flowchart  -  135 

C.7  "CORRCT"  Subroutine  Flowchart  -  136 

C. 8  "BON DRY"  Subroutine  Flowchart  -  139 

C. 9  "SRFLCT"  Subroutine  Flowchart  -  142 


A 


i 


j 

* 


A 


Speed  of  Sound 


*A  Denotes  the  value  of  a  variable  at  the 

node  to  the  left  of  a  discontinuity,  * 
can  be  any  variable  name 

AR  The  ratio  of  sound  speed  across  a  shock. 

A/B  (shock  moving  right) ,  B/A  (shock 
moving  left) 

*B  Denotes  the  value  of  a  variable  at  the 

node  to  the  right  of  a  discontinuity 

BDRY  3  denotes  left  boundary,  2  the  right 

boundary 

*BD  Variable  at  phantom  node 

CNTACT  3  digit  variable  denoting  contact  sur¬ 
face  location,  direction  of  travel,  and 
if  it  crosses  a  node 

COUNT  Counter  for  graphics  routines 

CSDIR  Contact  surface  direction,  2  to  the 

right,  3  to  the  left 

CSRMN  Riemann  variable  change  across  a  contact 
surface 

D2  Density  ratio  at  first  node  inside 

boundary 

DARRAY  Array  of  density  for  plotting 

DELAH  Change  in  A  from  I  to  1+1 

DELAL  Change  in  A  from  1-1  to  I 

DELQQH  Change  in  QQ  from  I  to  1+1 

DELQQL  Change  in  QQ  from  1-1  to  I 

Change  in  RR  from  I  to  1+1 


DELRRH 


DELRRL 
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DELT 

Time  step 
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DELTEX 
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DELTWL 

The  time  for  the  shock  to  reach  the  wall 
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Change  in  Q  from  I  to  1+1 
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Change  in  Q  from  1-1  to  I 

As 

DELX 

Interpolation  distance  (LMD+DELT) 

p 

DENS 

Density 

DLCD 

Density  to  the  left  of  the  contact 
discontinuity  in  the  exact  solution 

DLSH 

Density  to  the  left  of  the  shock  in  the 
exact  solution 

DLTA** 
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DQ 

The  jump  in  velocity  across  the  shock 
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calculation 
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The  jump  in  QQ  across  the  shock  calcu¬ 
lated  analytically  as  a  function  of  w 
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Gamma  (ratio  of  specific  heats) 
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For  graphical  output,  0  -  none  (tabular) , 
1  *  plots  all  variables,  2  *  compares 
density  with  exact  solution 

G1  1/ (G-l) 

G2  2/ (G-l) 
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HALT  Terminates  program  if  1,  set  by  condi¬ 

tions  not  coded 

INTEG(K)  Result  of  integrating  Z(K) 

**INT  Value  of  **  interpolated  between  nodes 
on  the  current  time  level 

12  Number  of  the  node  to  the  right  of  a 

discontinuity 

JSTOP  Number  of  time  levels  to  be  calculated 

LBDDR  Left  boundary  density  ratio 

LBDDRI  Left  boundary  density  ratio  at  time  zero 

LBDPR  Left  boundary  pressure  ratio 

LBDPRI '  Left  boundary  pressure  ratio  at  time  zero 
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left  boundary,  1  denotes  adjustable  pres¬ 
sure  at  the  left  boundary 

LBDTR  Left  boundary  temperature  ratio 

LBDTRI  Left  boundary  temperature  ratio  at  time 

zero 

LBNDRY  Denotes  left  boundary  condition,  open  or 
closed 

LMD(K)  The  characteristic  trajectories  in  the 
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LNODE  Array  of  left  most  node  to  be  corrected 

in  CORRCT 

Denotes  which  side  of  diaphragm  has  low 
pressure 
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The  measured  jump  in  QQ  across  the 
shock,  from  A  to  B 
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ND 

Double  precision  value  of  N 

NEW** (I) 

Stored  values  of  **  for  the  next 
level 
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PARRAY 

Array  of  pressures  for  plotting 
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The  ratio  of  the  pressure  across 
A/B  (right) ,  B/A  (left) 

a  shock. 

PRESS 
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PLTCNT 

*PRIM(K) 

P2 
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QARRAY 
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QQ 
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RBDPRI 


Pressure 

Initial  pressure  ratio  across  the  diaphragm 

Counter  for  graphics  routines 
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derivative  of  *  at  the  current  time  level 

Pressure  ratio  at  first  node  inside 
boundary 

Absolute  fluid  velocity 
Array  of  velocities  for  plotting 
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Initial  velocity  at  right  boundary 
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Initial  right  boundary  density  ratio 
Right  boundary  pressure  ratio 
Initial  right  boundary  pressure  ratio 
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RBDTR 

RBDTRI 

RBNDRY 

RNODE 

RR 

RXX 

S 

SAP 

SARRAY 

SAVG 

SBP 

SHKDIR 

SHOCK 

SIGMA 

SK 

SKIP 


Value  of  0  denotes  a  constant  pressure  at 
right  boundary,  while  1  is  for  adjustable 
pressure 

Right  boundary  temperature  ratio 

Initial  right  boundary  temperature  ratio 

Denotes  right  boundary  open  or  closed 

Array  for  right  most  node  to  be  corrected 
in  CORRCT 

Q-A*S  (extended  Riemann  variable) 

Node  defining  the  right  interval 
Entropy 

Entropy  to  the  left  of  the  shock  for 
flows  right  or  entropy  to  the  right  of 
the  C.S.  for  flows  left 

Array  of  entropy  for  plotting 

Average  entropy 

Entropy  to  the  right  of  the  C.S.  for 
flows  right  or  entropy  to  the  left  of 
the  shock  for  flows  left 

Shock  direction  of  travel,  3  to  the  left, 

2  to  the  right 

3  digit  variable  denoting  shock  location, 
direction  of  travel,  and  if  it  crosses  a 
node 

Spatial  location  of  discontinuities 
Sigma (L,J)  where  L  indicates  the  type  of 
discontinuity  and  J  indicates  the  time 
level:  1 — current  level,  2 — level  being 
calculated 

Integer  that  denotes  relative  location  of 
shock  near  boundaries 

Variable  which  indicates  how  many  time 
steps  between  calls  to  output  routines 


**STEP 


The  change  in  time  of  **  at  a  node  used 
to  step  up  to  the  next  time  level 


I 

T 

Time  since  initial  conditions 

si  T 

TEMP 

Temperature 

■ 

TRI 

Initial  temp,  ratio  across  the  diaphragm 

9 

TS 

Time  for  shock  to  travel  one  interval 

I 

T2 

Temperature  ratio  at  first  node  inside 
boundary 

R  UA 

UA 

Velocity  relative  to  the  shock,  left  side 

1  UB 

UB 

Velocity  relative  to  the  shock,  right 
side 

1 

VHEAD 

Velocity  of  the  head  of  the  expansion 
wave  for  the  exact  solution 

B 

VCDE 

Velocity  of  the  contact  discontinuity 
for  the  exact  solution 

1 

VS 

The  shock  speed  (positive  right, 
negative  left) 

i 

VSE 

Velocity  of  the  shock  for  the  exact 
solution 

1 

VTAIL 

Velocity  of  the  tail  of  the  expansion 
wave  for  the  exact  solution 

i  i  w 

W 

Mach  no.  relative  to  a  standard  shock 

H  w 

Vector  of  the  principle  variables,  Q,R,S 

H  x 

X 

Location  in  spacial  plane  (I-1)*H 

n 

XARRAY 

Array  of  spatial  positions  for  plotting 

I 

XEXACT 

Array  of  six  X  values  for  the  exact 
solution 

1 

XINIT 

Initial  position  of  discontinuity  for 
exact  solution  plotting 

i 

X2 

Location  of  node  to  right  of  discont. 
along  the  spacial  axis 

ii 

Y 

(N+l)/2 

if  ■ 

YEXACT 

Array  of  six  density  values  for  the 
exact  solution 

z 


Z(K) 


A 

6 


e,<j> 


Vector  of  the  right  hand  sides  of  the 
governing  equations 

Small  spatial  change 

Small  change  with  respect  to  time 

Flow  angles  with  respect  to  reference 
coordinate  planes 
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The  investigation  of  wave  rotors1  and  their  possible 
application  in  gas  turbine  engines  at  the  Naval  Postgraduate 
School's  (NPS)  Turbopropulsion  Laboratory  prompted  the 
development  of  a  one-dimensional  unsteady  flow  computer 
code.  The  QAZ1D  method  developed  by  Verhoff  [Ref.  1]  was 
chosen  over  other  methods  [Ref.  2]  for  reasons  which 
included  the  following: 

1.  The  method  is  based  on  the  use  of  characteristics. 
Such  methods  can  model  wave  propagation  accurately. 

2.  The  use  of  a  natural  streamline  coordinate  system 
eases  the  difficult  task  of  computing  with  two  and 
three  dimensional  grids. 

3.  The  equations  are  written  in  a  form  which  allows  a 
straightforward  extension  to  viscous  flows. 

Salacka  [Ref.  2]  verified  the  development  of  the  equations 

and  implemented  a  one-dimensional  Euler  solution  of  the 

i 

-  1The  wave  rotor  consists  of  simple  tube-like  passages 
along  and  around  the  periphery  of  a  drum  which  rotates 
between  end  walls  containing  inlet  and  outlet  (partial 
admission)  gas  ports.  Compression  and  expansion  processes 
are  arranged  to-  occur  in  a  cyclically  periodic  unsteady 
process  within  the  rotor  such  that  a  high  pressure  driver 
gas  is  used  to  compress  a  low  pressure  driven  gas.  The 
compressed  driven  gas  is  led  from  the  outlet  port  to  the 
combustor  and  brought  back  into  the  rotor  as  the  driver. 
The  expanded  driver  gas  exits  an  outlet  port.  Thus  the 
rotor  functions  as  both  turbine  and  compressor  with  direct 
gas-gas  energy  exchange  through  unsteady  wave  processes. 
The  technology  of  such  devices  requires  the  design  of 
appropriate  cycles  using  one-dimensional  calculations,  and 
analysis  of  the  multi-dimensional  unsteady  flows  such  as 
occur  during  port  opening  and  closing. 


shock  tube  problem  in  the  EULER1  code.  The  EULER1  code  was 
limited  to  high  pressure  set  on  the  left,  with  flow  and 
shock  velocities  to  the  right.  The  shock  wave  was  tracked 
and  corrected  for  in  the  flow,  but  the  contact  discontinuity 
was  not,  and  the  expansion  wave  was  not  tracked.  The  code 
also  did  not  treat  end  boundary  conditions. 

The  present  work  extended  the  EULER1  code  of  Salacka  to 
become  the  "Euler  ID  Version  2"  or  E1DV2  code  which: 

1)  can  handle  both  right  and  left  traveling  flow  and 
discontinuities 

2)  implements  tracking  and  correction  of  the  contact 
discontinuity 

3)  implements  tracking  of  the  expansion  wave 

4)  models  closed  wall  and  open  end  boundary  conditions 

5)  models  shock  and  contact  surfaces  within  a  grid 
interval. 

Section  II  and  Appendix  A  describe  the  analysis  underlying 
the  revisions.  The  Fortran  code,  E1DV2 ,  is  detailed  in 
Section  III  and  Appendix  C  and  graphical  results  of  four 
test  cases  to  demonstrate  the  new  features  of  the  code  are 


given  in  Section 

IV. 

A  discussion 

of  these 

results 

and 

limitations 

of 

the 

present  code  follow 

in 

Section 

V. 

Conclusions 

and 

recommendations  are 

made 

in 

Section 

VI. 

Appendix  B  contains  instructions  to  operate  the  code  on  the 
NPS  computer  and  the  FORTRAN  listing  is  included  as  Appendix 


II.  ANALYTICAL  AND  COMPUTATIONAL  APPROACH 


A.  QAZ1D  DESCRIPTION 

The  QAZ1D  analysis  and  numerical  solution  scheme  is  an 
explicit,  non-conservative  method  that  uses  extended  Riemann 
variables  to  model  the  multi-dimensional  Euler  equations  in 
a  natural  streamline  coordinate  system.  Table  I  lists  a 
summary  of  the  applicable  equations  which  are  derived  in 
detail  in  [Ref.  2]. 

1.  The  Coordinate  System 

Figure  2.1  shows  the  natural  streamline  coordinate 
system  (s,n,m)  relative  to  a  fixed  rectangular  cartesian 
system  (x,y,z).  The  right-hand  orthogonal  system  moves  in 
curvilinear  translation  along  a  streamline.  The  "s" 
direction  is  measured  in  the  direction  of  the  flow,  tangent 
to  the  streamline.  The  "n"  direction  is  perpendicular  to 
the  s  direction  in  a  plane  perpendicular  to  the  x-z  plane. 
The  "m"  direction  is  normal  to  the  plane  containing  the  n 
and  s  directions.  The  angle  p  is  the  angle  between  the  s-n 
and  x-y  planes,  and  the  angle  9  is  the  angle  between  the 
velocity  vector  and  the  x-z  plane  in  the  s-n  plane.  In  the 
one-dimensional  problem  treated  here,  the  s  direction  is 
such  that  velocity  to  the  right  is  positive,  and  to  the  left 
is  negative.  [Ref.  l:p.  2] 


TABLE  I 


SUMMARY  OF  THE  EULER  EQUATIONS 
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Figure  2.1  The  Natural  Coordinate  System 


2.  Extended  Riemann  Variables 

The  "Riemann  Variables”  in  most  publications  are 
defined  as 


Ri  "  <j  -  <7tt)A 

and  (2.1) 

R2  -  «  + 

where  q  is  the  velocity  magnitude,  A  is  the  speed  of  sound, 
and  y  is  the  ratio  of  the  specific  heat  at  constant  pressure 
to  the  specific  heat  at  constant  volume  for  a  particular  gas 


[Ref.  3:p.  5].  Verhoff  modified  the  definition  to  obtain 
"Extended  Riemann  Variables,"  defined  as 

Q  =  q  +  AS  (1.4) 

R  =  q  -  AS  (1.5) 


where  S  is  the  modified  entropy  given  in  Table  I  [Ref.  l:p. 

2]. 

3 .  The  Governing  Ecmations 

The  Euler  equations  in  Table  I  were  developed  in 
detail  [Ref.  2:Appendix  A]  from  the  equations  for 
conservation  of  mass,  momentum,  and  energy  for  -a  control 
volume.  The  underlying  assumptions  comprise  inviscid  flow, 
negligible  body  forces,  no  heat  transfer,  and  a  perfect  gas. 
These  assumptions  and  the  definition  of  modified  entropy 
were  used  in  transforming  the  equations  into  a  form  suitable 
for  solution  using  the  method  of  characteristics. 

B.  CHARACTERISTIC  SOLUTION 

Equations  (1.6),  (1.7),  (1.8),  (1.9)  and  (I. 10)  are  a 

system  of  slightly  coupled  quasi-one-dimensional  partial 
differential  equations  (PDE)  with  eigenvalues  q+A ,  q,  and  q- 
A.  These  equations  can  be  rewritten  along  the  characteris¬ 
tics  in  the  time-space  domain  to  become  ordinary 
differential  equations  (ODE)  which  are  solved  by  quadrature 
and  interpolation.  The  propagation  of  each  characteristic 


expressed  in  matrix  form  by  setting 


q+A  0 


0 


q-A  0 


-  ^A(s  -~j)  [q 


Then  equation  (2.2)  becomes,  for  the  equation  set, 


f +  1X1 1  •  z  • 


Figure  2 . 2  shows  the  characteristic  curve  in  the  space- 
time  domain  between  two  nodes  within  the  computational  mesh. 
Equation  (2.2)  has  a  solution 


t+dt  _ 

5w  =  -  Aw  +  /  z  dt 

t 


(2.5) 


where  fiw  is  the  change  due  to  time  at  a  fixed  location,  and 
Aw  is  the  change  due  to  displacement  at  a  fixed  location, 
along  the  characteristic  curve.  The  value  of  Aw  is 
calculated  by  interpolation  through  an  iterative  scheme 
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Figure  2.2  Characteristic  Curve  Within  Computational 
Mesh 


using  initial  guesses  for  the  characteristic  slope,  X ,  from 
values  at  time  level  t.  The  line  integral  is  approximated 


t+dt 


B  r 

z  dt  *  /  y  ds 
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/  'fdS 


f(Si-(Si-AS)) 
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The  spatial  derivatives  of  q  and  A  in  T  are  estimated  from 
values  at  A  and  s^.  Finally,  9w  is  calculated,  and  each 
node  is  updated  to  the  next  time  level. 

The  Courant-Friedrichs-Lewy  (CFL)  convergence  condition 
requires  the  domain  of  dependence  of  the  numerical  approxi¬ 
mation  to  include  the  domain  of  dependence  of  the  differen¬ 
tial  equation  [Ref.  4:pp.  336-3371.  Figure  2.3  shows  that 


the  computational  grid  for  the  QAZ1D  method  satisfies  this 
condition  by  having  interpolated  initial  data  points  (points 
1,  2,  and  3)  in  between  previous  solution  nodes  (points  4, 
5,  6)  .  The  q+A  characteristic  is  estimated  from  values  at 
point  4,  q  characteristic  from  values  at  point  5,  and  q-A 
characteristic  from  values  at  point  6.  The  use  -of  a  maximum 


time  step  keeps  the  domain  of  dependence  of  the  numerical 
scheme  just  outside  that  of  the  actual  equations,  aiding 
convergence.  Unlike  finite  difference  schemes  that  are 
conservative  and  require  a  numerical  diffusion  or  viscosity 
to  handle  oscillations  or  smearing,  this  method  is  stable 
without  numerical  smoothing  or  artificial  dissipation  [Ref. 
5:pp.  562-583].  However,  the  equations  do  not  account  for 
the  discontinuities  across  shocks  and  contact  surfaces 
properly.  Thus  after  each  time  step  a  one-point  correction 
technique  is  used  to  correct  for  shocks  and  contact 
surfaces. 

C.  DISCONTINUITIES 

1.  Exoansion/Rarefaction  Waves 

The  head  and  tail  of  the  expansion  wave  travel  along 
the  characteristic  as  gradient  discontinuities.  Flow 
velocity,  speed  of  sound,  pressure,  and  density  are  continu¬ 
ous  across  a  gradient  discontinuity  but  the  spatial  deriva¬ 
tives  are  discontinuous.  Entropy  remains  constant  through 
the  rarefaction  wave.  The  head  of  an  expansion  wave  moves 
into  undisturbed  flow  with  a  speed  q+A,  while  the  tail  of 
the  expansion  wave  has  a  velocity  q-A.  No  special  treat¬ 
ment  of  these  waves  is  required;  the  characteristics  method 
accurately  tracks  and  models  their  influence.  Expansion 


waves  are  generated  when  two  shocks  collide,  a  diaphragm  is 
burst,  a  contact  surface  and  shock  intersect,  a  shock  leaves 


an  open  boundary,  or  an  open  boundary  has  a  pressure  lower 
than  the  pressure  inside  the  tube.  [Ref.  3:pp.  13-15] 

2.  Shocks 

Shocks  are  created  when  a  diaphragm  is  burst, 
compression  waves  generated  at  an  open  boundary  coalesce 
into  a  wavefront  with  sufficient  strength,  or  when  a  shock 
and  contact  surface  collide.  Pressure,  velocity,  density, 
and  entropy  are  discontinuous  across  a  shock  [Ref.  6:pp.  so¬ 
il].  The  equations  of  motion  (I. 6) -(I. 10)  are  not  correct 
for  calculating  changes  through  shocks.  A  method  to  correct 
for  a  shock  is  to  use  the  ratio  of  the  jump  in  the  extended 
Riemann  variable  across  the  shock  to  the  speed  of  sound 
downstream  of  the  shock  to  find  the  incoming  Mach  Number 
relative  to  the  shock  (w)  [Ref.  l:p.  3].  Figure  2.4  shows 
the  case  of  a  shock  headed  to  the  right.  The  shock  velocity 
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Figure  2.4  Shock  Wave  with  High  Pressure  on  Left 


(Vs)  is  imposed  on  the  flow  to  produce  a  stationary  shock 
condition.  Flow  variables  are  then  corrected  using  normal 
shock  relations.  Figure  2.4  applies  to  a  condition  where 
high  pressure  is  on  the  left.  The  shock  is  moving  to  the 
right  at  Vs  into  flow  with  velocity  qB.  Table  II  lists  the 
equations  developed  by  Salacka  [Ref.  2]  for  relating  the 
extended  Riemann  variable  change  to  Mach  Number.  Velocities 
are  non-dimensional ized  by  the  speed  of  sound  downstream  of 
the  shock  direction,  AB,  and  pressures  and  densities  by  the 
respective  values  downstream. 

Figure  2.5  depicts  equation  (II. 2)  over  a  range  of 
Mach  Numbers  from  1.0  to  4.0.  The  curve  is  approximated  by 
the  quadratic  equation 


qa~9b  = 
ab 


-2.7574  4-  3.1573W  -  0.2863W2 


(2.7) 


If  the  high  pressure  is  to  the  right  the  shock  moves 
as  Fig.  2.6  illustrates,  and  the  governing  equations  are  as 
listed  in  Table  III.  Again,  downstream  conditions  are  used 
to  non-dimensionalize  velocity,  pressure,  density,  and  the 
Riemann  variable  change.  The  development  of  the  equations 
in  Table  III  is  given  in  Appendix  A.  A  plot  of  equation 
(III.  2)  over  a  range  of  w  from  1.0  to  4.0  is  shown  in  Fig. 
2.7.  The  curve  is  approximated  by  the  quadratic  equation 
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=  2.7574  -  3.15732W  +  0.2863W2 


(2.8) 


TABLE  II 

SHOCK  EQUATIONS  FOR  HIGH  PRESSURE  ON  LEFT 


Definitions 

Relative  Incoming 
Mach  Number:  w 


AB 


— V 


Ecruations 

Q  Riemann 
Variable  Change: 
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TABLE  II  (CONTINUED) 


Velocity  Change: 
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Entropy  Change: 
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Figure  2.6  Shock  Wave  with  High  Pressure  on  Right 


TABLE  III 

SHOCK  EQUATIONS  FOR  HIGH  PRESSURE  ON  RIGHT 


Relative  Incoming 
Mach  Number: 


q.  -V 
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(IH.l) 


Equations 

R  Riemann 
Variable  Change: 
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TABLE  III  (CONTINUED) 
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Figure  2.7  Extended  Riemann  Variable  Change  with  Mach  Number,  High 
Pressure  on  Left 


An  examination  of  equations  (II. 2)  and  (III. 2) 


reveals  that 


V*A  fQA"% 

—£ - " 


i 


i 


This  is  expected  since  in  flow  to  the  right,  the  Q  extended 
Rienann  variable  is  associated  with  the  q+A  characteristic. 
Also,  since  velocity  is  defined  as  positive  to  the  right,  q 
is  positive.  But  for  flow  to  the  left,  velocity  is  defined 
as  negative,  thus  q  is  negative.  Moretti  [Ref.  3]  showed 
that  the  downstream  running  Riemann  variable 
(q+A  characteristic)  has  a  lower  magnitude  of  discontinuity 
across  a  normal  shock  than  the  upstream  running  Riemann 
variable  (q-A  characteris-tic) ,  and  is  thus  preferable.  In 
flow  to  the  left,  the  Riemann  variable  associated  with 
downstream  is  the  R  Riemann  variable,  though  the  character¬ 
istic  associated  with  the  R  Riemann  variable  is  q-A ,  in  flow 
to  the  left  q  is  negative  and  thus 


Therefore,  in  fluid  traveling  left  the  R  Riemann  variable  is 
used  to  find  Mach  Number,  and  in  fluid  traveling  right  the  Q 
Riemann  variable  is  used. 


If  the  shock  is  between  two  nodes  with  no  contact 
surface  within  the  same  interval,  the  method  to  calculate 
the  Mach  Number,  w,  is  as  follows: 

1)  The  change  in  the  appropriate  extended  Riemann  varia¬ 
ble  across  the  interval  is  measured;  R  variable  for 
flow  left,  and  Q  variable  for  flow  right.  Call  this 
change,  Am. 

2)  Use  Am  in  equation  (2.7)  or  equation  (2.8)  for  the 
appropriate  flow  to  calculate  w. 

3)  With  equation  (II. 2)  or  equation  (II. 3),  for  flow 

right  or  left  respectively,  use  w  to  calculate  the 

analytical  Riemann  variable  change,  Ae. 

4)  If  this  exact  value  of  Ae  is  not  equal  to  Am  then 

calculate  a  new  guess,  Ag  from 

Agi+1  =  Ag-J_  +  (Am  -  Ae)  (2.10) 

and  then  repeat  steps  2  to  4  until  Ae  equals  Am. 

Thus  the  correct  value  of  w  is  determined,  and  shock 
corrections  from  the  normal  shock  relations  are  valid.  Vs 
is  determined  from  equation  (II. 1)  or  equation  (III.l) 

depending  on  flow  direction.  Flow  direction  will  be 
initially  determined  by  the  side  with  high  pressure.  High 
pressure  on  the  left  means  flow  will  travel  to  the  right. 
Likewise,  high  pressure  on  the  right  means  flow  will  travel 
to  the  left. 

3 .  Contact  Surface 

A  contact  surface  is  a  boundary  between  regions  of 
flow  which  are  different  in  composition  or  which  have 
undergone  different  thermodynamic  histories.  Thus  a  contact 
surface  is  formed  when  the  diaphragm  is  burst  in  a  shock 


tube  separating  the  gas  initially  in  the  driver  from  the  gas 


initially  in  the  driven  tube  [Ref.  6:pp.  23-29].  In  a  wave 


rotor,  a  contact  surface  would  separate  the  combustor  gases 


from  the  cooler  inlet  air.  Contact  surfaces  are  also  caused 


by  two  shocks  of  opposite  families  colliding  or  a  shock  over 


taking  a  shock  of  the  same  family  [Ref.  3:pp.  15-18]. 


Figure  2.8  illustrates  the  changes  in  physical 


properties  through  the  contact  surface  in  the  shock  tube 


problem.  The  gas  to  the  right  has  been  compressed  and 


heated  by  the  shock  wave,  with  that  to  the  left  cooled  and 
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Figure  2.8  Behavior  of  Physical  Properties  Through 
a  Contact 'Surface 


expanded.  The  density,  speed  of  sound,  and  modified  entropy 
are  discontinuous  across  the  contact  surface. 


Since  pressure  is  constant  across  the  contact 
surface,  using  the  perfect  gas  equation  of  state  and  the 
first  law  of  thermodynamics  it  is  shown  in  Appendix  A  for 
the  shock  tube  problem  that 
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(2.11) 


Since  the  flow  behind  the  contact  surface  to  the 
tail  of  the  rarefaction  wave  is  isentropic,  the  value  of  SB 
will  be  known.  The  value  of  is  the  value  of  entropy 
behind  the  shock.  Subtracting  equations  (1.4)  and  (1.5) 
from  each  other  and  dividing  by  entropy  SB,  A  B  can  be 
determined.  Thus  the  unknown  speed  of  sound  A a  is  obtained 
from  equation  (2.11).  With  qB  equal  to  qA,  and  SA  the 
Riemann  variables  just  behind  the  contact  surface  are  calcu¬ 
lated  from  equations  (1.4)  and  (1.5).  For  a  contact  surface 
traveling  to  the  left  it  is  only  necessary  to  interchange 
the  subscripts.  The  velocity  of  the  contact  surface  is 
easily  determined  since  it  moves  at  velocity  q  along  the  q 
characteristic  line  associated  with  S.  [Ref.  3:p.  16] 

4 .  Contact  Surface/Shock  Interaction 

When  the  shock  and  contact  surface  are  within  the 
same  interval  the  calculation  of  the  shock  incoming  Mach 
Number,  w,  must  be  modified.  Figure  2.9  shows  a  plot  of  the 


Figure  2.9  Modified  Entropy  Distribution  for  Shock 
and  Contact  Surface  Traveling  Right 

modified  entropy  change  across  the  interval.  The  use  of  the 
Riemann  values  art  the  nodes  to  estimate  the  Riemann  variable 
change  across  the  shock  would  be  incorrect.  The  correct 
Riemann  variable  change  would  be,  for  flow  to  the  right, 


ASK  = 


(2.12) 


The  correct  Mach  Number,  w,  is  determined  using  the 
technique  of  Moretti  modified  for  QAZ1D  [Ref.  3:pp.  23-24]. 

Referring  to  Fig.  2.9,  the  technique  for  flow  to  the 
right  is  as  follows: 

1)  Calculate  w  as  if  no  contact  surface  were  in  the 
interval.  With  values  at  node  i+1  known,  determine 
S3.  Then  use  w  and  S^  in  equation  (II. 7)  to  solve  for 
S^.  This  is  the  initial  estimate 'of  Ss.  Let  SB> =  Ss. 

2)  The  Q  Riemann  Variable  change  across  a  contact  surface 
is  shown  in  Appendix  A  to  be  given  by 
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3)  Define 


qa'-Ob 


(2.14) 


then 


Am  ■  ASK  +  ACS 


(2.15) 


can  be  solved  to  obtain  ASK. 

4)  Using  ASK  as  Am  in  the  shock  iteration  scheme  given  in 
Section  II. C. 2  calculate  a  new  w. 
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5)  With  this  new  w  and  Sq  use  equation  (II. 7)  to  obtain  a 
new  SA. 

6)  Compare  SA  with  Sg.  If  they  are  not  egual  to  within 
an  acceptable  error,  calculate  a  new  estimate  for  Sq 
using 
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(2 . 16) 
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Iterate  until  convergence  is  achieved. 

This  will  result  in  the  proper  w  for  a  condition 
with  shock  and  contact  surface  interacting. 

For  flow  to  the  left,  it  is  only  necessary  to  inter¬ 
change  the  subscripts  A  and  B,  and  use 
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for  the  extended  Riemann  Variable  change, 
illustrates  the  left  traveling  condition. 


Figure  2.10 


Figure  2.10  Modified  Entropy  Distribution  for  Shock 
and  Contact  Surface  Traveling  Left 

D.  BOUNDARY  CONDITIONS 
1.  Open  Boundary 

At  the  open  boundary,  a  reference  pressure  ratio, 
Pref  *  P«/PA  specified,  where  P*  is  the  pressure  to  be 
held  at  the  boundary,  and  PA  is  the  pressure  just  inside  the 
tube.  Only  the  case  where  P*,  <  PA  is  considered.  This 
means  that  an  outflow  condition  at  the  boundary  will  result 
when  Pffl  <  PA,  with  an  expansion  wave  traveling  in  from  the 
boundary.  Thus  locally  at  the  boundary,  conditions  are 
isentropic. 

The  computational  grid  for  the  left  boundary  is 
shown  in  Fig.  2.11.  A  phantom  node,  L,  is  used  outside  the 
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Figure  2.11  Left  Boundary  Computational  Grid 


computational  mesh  to  enable  simple  enforcement  of  constant 


pressure  and  entropy  at  the  mesh  boundary  node,  l. 


The  approach  is  to  determine  the  required  variables 


at  node  L,  then  transfer  these  values  onto  node  1,  and  vary 


q  at  the  boundary  to  meet  the  required  boundary  conditions. 


First,  using  Pref  at  node  L  and  S  at  node  2  in  equation 


(1.3)  rewritten  in  the  form 


(S  --ir)  (-y(y-D)  -1/y 
Or  "  ( (1/Pref )  (exp  )) 


(2.18) 


Or,  the  density  ratio  at  the  boundary,  is  calculated.  Using 
the  perfect  gas  equation  of  state  in  equation  (2.18)  the 
temperature  ratio,  Tr  is  given  by 


=  (o 


(Y-l) 

R 


)  (exp 


(Y(1-Y)(S 


) 


(2.19) 


and  can  be  calculated  once  Pr  is  known.  With  an  initial 
velocity  at  the  boundary,  q^,  assumed  the  Riemann  variables 
R  and  Q  are  determined  for  node  L  using 


R  = 


(2.20) 


and 


Q  -  qb  +  ^"Tr  S  (2.21) 

Subtracting  equation  (2.20)  from  (2.21),  and 
dividing  by  twice  S  yields  A  at  node  L,  where  A  =  /y  RGT. 
The  values  of  R,  Q,  A,  q^,  and  S  are  now  known  at  node  L. 
These  values  are  assigned  to  node  1,  and  the  QAZ1D  algorithm 
is  applied  to  the  boundary  node  as  if  it  were  an  interior 
node.  The  result  is  an  estimate  of  the  new  values  for  R,  Q, 
and  S  at  node  1.  These  are  used  to  obtain  a  new  estimate  of 
the  pressure  ratio  at  node  1,  call  it  Prn*  Prjj  is  compared 
with  Pref,  and  the  new  value  of  S  with  the  old  value  of  S  at 
node  1.  A  new  q^  is  calculated  by  solving  equations  (2.20) 
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and  (2.21)  simultaneously.  The  process  is  repeated  until 
entropy  and  pressure  remain  constant  at  the  boundary.  For 
the  right-hand  boundary,  Fig.  2.12  shows  the  computational 


•t 

right  grid  boundary 
B 


Figure  2.12  Right  Boundary  Computational  Grid 


grid  with  the  phantom  node,  r,  to  the  right  of  the  mesh 
boundary  node,  n.  The  procedure  for  the  right  boundary  is 
the  same  as  for  the  left  boundary. 

When  a  shock  travels  across  the  open  boundary,  the 
boundary  is  first  corrected  using  the  method  described  in 
Section  II. C. 2.  Subsequently,  the  pressure  is  returned  to 
give  the  constant  value  specified  for  Pref,  using  the  above 


procedure.  The  result  will  be  an  expansion  wave  traveling 
inward . 


For  supersonic  flow,  all  of  the  variables  are  deter¬ 
mined  from  values  at  the  interior  nodes  [Ref.  l:p.  3]. 

2 .  Closed  Boundary 

The  solid  boundary  is  imposed  by  setting  the 
velocity  at  the  boundary,  q,  equal  to  zero.  The  same 
computational  grid  is  used  as  that  for  an  open  boundary,  but 
a  different  technique  is  used  in  assigning  values  to  the 
phantom  node.  Referring  to  Fig.  2.11,  the  velocity  at  node 
2,  q2,  is  known.  By  setting 

SL  *  “  <22  (2.22) 

the  wave  impacting  the  boundary  will  be  met  by  a  wave  of 
equal  velocity  but  opposite  in  direction  and  will  result  in 
node  1  appearing  as  a  solid  boundary.  Further,  to 

facilitate  the  QAZ1D  method,  set 


This  will  enable  the  computation  of  the  new  boundary  node  1 
values  as  an  interior  point  with 


«» 
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qi  -  o 

The  right  boundary  is  handled  in  the  same  way. 

Figure  2.13  illustrates  the  sequence  of  events  when 
a  shock  reflects  off  a  solid  boundary  on  the  right  [Ref.  7: 


A)  Shock  Approaches  Solid  B)  Shock  Reflects  Off 

Boundary  Solid  Boundary 

Figure  2.13  Shock  Reflecting  at  Solid  Boundary 

pp.  172-195].  Since  the  velocity  behind  the  reflected 
shock,  qB,  is  zero,  and  the  velocity  in  front  of  the  shock, 
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Aq  =  (qs  -  3a)/AA 


(2.27) 


Rearranging  equation  (III. 6),  and  using  Y=  1.4  it 


can  be  shown  that 


w  =  Jl  +  0 . 3 6 ( Aq) 2  -  0 . 6Aq 


(2.28) 


This  is  an  exact  analytical  value  for  w  based  on  stationary 
normal  shock  relations.  The  speed  of  sound  ratio,  Ag/AA/  is 
calculated  from  equation  (III. 3)  and 


ab  -  Aa(Ab/Aa) 


(2.29) 


Entropy  behind  the  shock,  Sq,  is  determined  from  equation 
(III. 7)  since  is  known.  Then  the  extended  Riemann 

variables,  R  and  Q,  behind  the  reflected  shock  are 
calculated  from  equations  (1.4)  and  (1.5)  respectively. 

For  a  shock  reflected  from  the  left  boundary,  the 
subscripts  A  and  B  are  interchanged  and  using  equation 
(II. 6)  the  Mach  Number  is  given  by 


w  =  /l  +  0.36(£q)2  +  0.6A' 


(2.30) 


The  equations  from  Table  II  are  used  instead  of 
those  from  Table  III  (used  in  the  case  of  reflection  from 
the  right  boundary) ,  since  the  shock  is  now  traveling  from 
left  to  right. 


III.  FORTRAN  PROGRAM  E1DV2 

A.  GENERAL  DESCRIPTION 

Prograir  E1DV2  is  a  Fortran  computer  program  which 
calculates  1-dimensional  unsteady  flow  based  on  the  QAZ1D 
formulation  and  method  of  solution  of  the  Euler  equations  of 
motion.  The  program  has  provisions  for  tracking  in  time  the 
location  of  shock  waves,  contact  surfaces  and  head  and  tail 
of  the  expansion  waves.  Its  intended  application  to  wave- 
rotc*r  design  and  analysis  explains  the  present  definition  of 
node  points  and  boundary  conditions.  In  its  current  form, 
the  code  describes  flow  in  a  tube  initiated  by  a  diaphragm 
bursting.  The  location  of  the  diaphragm,  whether  the  flow 
is  to  the  left  or  the  right,  and  whether  the  ends  of  the 
tube  are  closed  or  open,  can  be  varied. 

The  E1DV2  program  is  written  in  FORTRAN  77,  and  was 
developed  using  the  IBM  3033  System  370  computer.  The  use 
of  DISSPLA  subroutines  in  plotting  results  requires 
compiling  and  executing  using  VS  FORTRAN  procedures. 
Appendix  B  describes  the  procedures  for  running  the  E1DV2 
program  on  the  Naval  Postgraduate  School  IBM  computer 
system.  A  listing  of  the  code,  which  contains  its  own 


description  in  comment  statements,  is  given  in  Appendix  0. 


1)  first  order  accuracy 

2)  double  precision  numerics,  except  for  single  precision 
in  graphical  output 

3)  linear  interpolation  and  extrapolation  algorithms 

4)  quadratic  polynomial  approximations  for  curve  fitting 

5)  non-dimensional  variable  input  and  output 

6)  structured  subroutine  format. 

Table  IV  lists  all  the  subroutines  called  from  the  main 
program  and  brief  descriptions  of  their  purposes.  The 
extensive  use  of  subroutines  was  intended  to  allow  modifica¬ 
tions  and  extensions  of  the  code  without  major 
restructuring . 

B.  THE  MAIN  PROGRAM 

The  main  program  flowchart  is  illustrated  in  Appendix  C, 
Fig.  C.l.  The  conventions  and  a  comprehensive  list  of 
variables  are  listed  in  the  beginning  of  the  program.  The 
number  of  grid  nodes  is  set  by  the  user,  and  must  be  an  odd 
number.  Arrays  are  dimensioned  to  the  number  of  nodes  in 
the  main  program  prior  to  compiling  and  executing. 

Initial  values  for  temperature,  pressure,  and  density 
ratios  at  the  diaphragm  and  boundaries  are  entered  by 
editing.  Also,  the  initial  velocity  distribution  on  each 
side  of  the  diaphragm  is  set,  and  as  well  as  the  ratio  of 
specific  heats.  The  side  with  the  high  pressure  is  identi¬ 
fied.  Boundaries  are  set  either  closed  or  open.  In  the 


TABLE  IV 


TIME 

TRAK 

SWEEP 

CONDI 

COND2 

COND3 

COND4 

COND5 

COND6 

COND7 


SUBROUTINES 


Calculates  minimum  time  step 

Tracks  discontinuities,  calculates  shock  mach 
number,  w 

Advances  node  values  to  next  time  step 

Calculates  interpolated  values  of  Riemann 
variables,  3q/3s,  and  3A/ 3s  when  no  shock 
or  contact  surface  within  an  interval 

Calculates  interpolated  values  of  Riemann 
variable,  3q/3s,  and  3A/3 s  when  discontinuity 
within  interval  to  right  of  node  or  flow  is 
supersonic  to  the  right 

Calculates  interpolated  values  of  Riemann 
variable,  3q/3s,  and  3A/3s  when  discontinuity 
within  interval  to  left  of  node  or  flow  is 
supersonic  to  the  left 

Called  when  node  is  jumped  by  discontinuity 
to  load  node  information  into  array  for  use 
in  subroutine  CORRCT 

Calculates  interpolated  values  of  Riemann 
variable,  3q/3s,  and  3A/3s  when  discontinuity 
on  either  side  of  node  and  flow  is 
supersonic 

Would  contain  subroutine  to  advance  node 
when  contact  surface  is  in  front  of  shock 
after  intersecting  and  traveling  in  same 
direction.  Not  currently  in  use 

Would  contain  subroutine  to  solve  precisely 
when  and  where  shock  and  contact  surface 
would  intersect  within  an  interval.  Not 
currently  in  use 

Would  contain  subroutine  to  solve  precisely 
when  and  where  shock  and  contact  surface  would 
intersect  if  either  were  jumping  a  node  simul¬ 
taneously.  Not  currently  in  use 
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TABLE  IV  (CONTINUED) 


COND7S 

COND8 

CORRCT 

DELTAX 

SKJUMP 

CSJUMP 

EXTRAP 

INTERP 

DBURST 

BONDRY 

SP.FLCT 

BBDRY 

BORDER 

PLOT 

EXACT 


Would  contain  the  subroutine  to  solve  the 
Riemann  problem  of  shock  and  contact  surface 
intersecting.  Not  currently  in  use 

Would  contain  subroutine  to  advance  node  after 
shock  and  contact  surface  have  crossed  and  are 
moving  apart.  Not  currently  in  use 

Advances  nodes  jumped  by  discontinuity  to  next 
time  step 

Determines  relative  location  of  discontinuity 
within  an  interval 

Calculates  variable  change  across  a  normal 
stationary  shock 

Calculates  variable  change  across  a  contact 
surface 

Extrapolates  entered  values 

Interpolates  entered  values 

Calculates  initial  Riemann  values  at  node 
where  diaphragm  is  located 

Calculates  values  to  update  left  and  right 
boundary 

Calculates  new  shock  parameters  when  shock 
reflects  off  solid  boundary 

Alternate  interpolating  scheme  near  boundary 

Sets  graphical  borders 

Plots  q,  S,  P,  and  p 

Plots  exact  p  versus  calculated  p 

Lists  calculated  values  in  files  8,  9,  10 


LIST 


open  case,  either  constant  pressure,  or  an  adjustable 
pressure  is  specified.  The  latter  case  in  effect  extends 
the  tube,  and  allows  discontinuities  to  disappear  out  of  the 
tube  without  creating  expansion  or  compression  waves.  Exact 
values  for  velocities  of  expansion  wave,  shock,  and  contact 
surface,  and  density  ratio  are  specified  by  the  user  in  the 
input  section. 

Output  is  controlled  by  setting  the  variable  GRAPHS  to 
either  0,  1,  or  2 .  A  zero  creates  three  files  which  contain 
data  calculated  at  the  time  the  call  to  subroutine  LIST  is 
made.  A  value  of  1  causes  plots  of  pressure,  density, 
velocity,  and  entropy  distributions  to  be  created.  When 
GRAPHS  equals  2  a  plot  of  the  exact  density  distribution 
along  with  the  calculated  density  distribution  is  made. 
Calls  to  output  routines  are  determined  by  the  parameter, 
SKIP,  the  number  of  time  steps  between  calls,  which  is 
specified  by  the  user. 

Program  termination  can  occur  for  several  reasons. 
First,  the  program  terminates  when  a  maximum  number  of  time 
steps,  JSTOP,  is  reached.  JSTOP  is  specified  by  the  user. 
Second,  the  program  terminates  if,  during  execution,  any  of 
the  following  conditions  are  met: 

1)  The  shock  intersects  with  the  contact  surface 

2)  The  pressure  at  any  open  boundary  is  higher  than  the 
pressure  at  the  first  node  inside  the  boundary 

2)  The  shock  exits  an  open  boundary. 


In  these  cases  a  message  is  displayed  on  the  terminal 


describing  the  reason  for  termination. 

C.  THE  SUBROUTINE  PROGRAMS 

1.  The  "DBURST"  Routine 

Figure  C.2  in  Appendix  C  shows  the  flowchart  for 
"DBURST".  The  assumption  that  a  shock  and  contact  surface 
are  formed  immediately  causes  oscillating  numerical  solu¬ 
tions  that  dampen  within  five  to  ten  time  steps.  By 
mathematically  "bursting"  the  diaphragm  prior  to  time  zero, 
a  solution  at  the  node  where  the  diaphragm  is  located  can  be 
determined  at  time  zero.  The  result  is  a  more  accurate 
representation  of  the  wave  structure  and  a  dampening  of  a 
reduced  transient  solution  within  three  time  steps. 
"DBURST"  calculates  the  Riemann  variables  that  would  exist 
behind  a  shock  and  contact  surface  that  have  moved  an 
infinitesimal  amount  after  the  diaphragm  burst,  and  assigns 
them  to  the  node  where  the  diaphragm  was  situated. 

2.  The  "TIME"  Routine 

The  "TIME"  subroutine  was  not  changed  from  that 
developed  by  Salacka  [Ref.  2:p.  35].  The  purpose  is  to 

compute  the  maximum  allowable  time  step  but  ensuring  that 
all  characteristic  trajectories  over  the  time  step  are 
within  one  interval  between  nodes.  The  time  step  is 
determined  from 
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DELT 


H/ABS  [  Q+A  ] 


(3.1) 


at  every  node,  and  the  minimum  time  step  is  used. 
3.  The  "TRAK"  Routine 


This  routine  computes  the  new  locations  of  the 
shock,  contact  surface,  and  head  and  tail  of  the  expansion 
wave  after  the  time  step,  DELT,  calculated  in  "TIME.”  In 
Appendix  C,  Fig.  C.3  shows  the  flowchart  for  the  "TRAK" 
subroutine.  The  shock  velocity,  Vs,  incoming  Mach  Number, 
w,  and  pressure,  density,  and  sonic  velocity  ratios  PR,  DR, 
and  AR  respectively  are  determined.  The  equations  for  flow 
left  and  right  are  the  same,  except  for  DQ,  the  velocity 
gradient,  which  is  different  only  by  sign.  Thus  by  proper 
bookkeeping  the  same  equations  are  coded  once  and  used  in 
either  direction.  The  initial  direction  of  the  shock  and 
contact  surface  are  set  from  code  that  interprets  which  side 
has  the  high  pressure  at  time  zero.  The  contact  surface 
speed  is  determined  after  the  shock  speed  is  established. 
The  situation  at  time  zero  is  handled  as  in  "DBURST"  to 
ensure  that  the  shock/contact  surface  interaction  is 
properly  modeled.  Since  there  is  a  certain  transient 
solution  that  decays  after  three  time  steps,  tracking  of  the 
expansion  wave  is  not  initiated  until  then. 

The  head  of  the  expansion  wave  travels  at  a  velocity 
corresponding  to  q+A  ,  and  the  tail  with  speed  q-A.  In  the 
case  of  flow  traveling  to  the  left,  this  is  reversed. 


"BONDRY"  subroutine  is  called  and  returns  the  updated  values 
at  those  nodes.  In  between,  the  correct  and  most  effective 
algorithm  to  update  the  node  is  determined.  These 
algorithms  are  coded  into  eight  different  condition 
subroutines.  These  routines  calculate  the  interpolated 
Riemann  variables  associated  with  the  three  characteristics 
(q+A,  q,  and  q-A  )  ,  plus  the  spatial  derivatives  at  3q/";S  and 
34/ 3s.  This  allows  "SWEEP”  to  calculate  Aw  and  2.  The 
integral  of  2  is  then  determined,  and  the  update  for  the 
variables  of  that  node  in  5w  are  computed.  3w  is  stored 
until  the  entire  mesh  has  been  swept.  Then  the  variable 
arrays  (w)  are  updated  to  the  next  time  interval. 

To  choose  the  proper  condition  routine,  the 
following  information  about  conditions  in  the  interval  on 
either  side  of  the  node  are  expressed  in  two  3  digit  integer 
variables,  SHOCK  and  CNTACT .  The  value  of  these  variables 
provide  a  code  for  the  following  information: 

SHOCK — The  existence  of  a  shock  within  an  interval,  and, 
if  so,  the  direction  of  travel,  and  if  the  shock 
will  cross  the  node  in  front  of  it  within  the 
current  time  step. 

CNTACT-The  existence  of  a  contact  surface  within  an 
interval,  and,  if  so,  the  direction  of  travel,  and 
if  the  contact  surface  will  cross  the  node  in  front 
of  it  within  the  current  time  step. 

Figure  3.2  illustrates  the  regions  around  a  node.  Table  V 

lists  the  values  that  SHOCK  and  CNTACT  may  have,  and  their 

meaning.  Figures  3.3,  3.4,  and  3.5  illustrate  examples  of 

what  different  SHOCK  and  CNTACT  values  mean.  At  each  node 
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Figure  3.2  Interval  and  Node  Description  used  in  "SWEEP 
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Figure  3.3  Schematic  of  SHOCK  =  331  and  CNTACT  =  232 


TABLE  V 


VALUES  OF  PARAMETERS  SHOCK  AND  CNTACT 


SHOCK 

Located  Direction  Crosses  node  in 


k&li&l 

of  travel 

100 

No  Shock 

N/A 

N/A 

321 

Left 

Right 

Yes 

322 

Left 

Right 

No 

331 

Left 

Left 

Yes 

332 

Left 

Left 

No 

221 

Right 

Right 

Yes 

222 

Right 

Right 

No 

231 

Right 

Left 

Yes 

232 

Right 

Left 

No 

CNTACT 

100 

No  contact  surface 

N/A 

N/A 

321 

Left 

Right 

Yes 

322 

Left 

Right 

No 

331 

Left 

Left 

Yes 

332 

Left 

Left 

No 

221 

Right 

Right 

Yes 

222 

Right 

Right 

No 

231 

Right 

Left 

Yes 

232 

Right 

Left 

No 

the  code  examines  the  shock  and  contact  surface  condition, 
along  with  relative  location  of  the  shock  to  a  contact 
surface  if  both  exist  near  the  node,  if  the  flow  is  super¬ 
sonic  or  subsonic  at  the  node,  and  the  direction  the  flow  is 
traveling.  The  appropriate  algorithm  is  determined  and 
called  in  the  form  of  a  condition  subroutine. 

The  "SWEEP”  routine  is  designed  to  handle  all 
possible  shock  and  contact  surface  combinations.  Because 
the  code  does  not  solve  the  situation  where  the  shock  and 
contact  surface  intersect,  as  illustrated  in  Fig.  3.6,  then 
any  combinations  after  the  shock  and  contact  surface  have 
crossed  call  condition  routines  that  currently  contain  only 
messages.  The  code  was  intentionally  prepared  to  be  easily 
extended  to  treat  the  intersection  of  the  shock  and  the 
contact  surface,  and  beyond. 
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Figure  3.6  Schematic  of  SHOCK  =  231  and  CNTACT 


The  ten  subroutines  that  are  called  by  ’•SWEEP"  to 
calculate  the  interpolated  values  of  the  extended  Riemann 
variables,  R  and  Q,  plus  S,  9q/3s,  and  ?A/3  s  are  CONDI, 
C0ND2 ,  C0ND3 ,  C0ND4 ,  C0ND5 ,  C0ND6 ,  C0ND7 ,  C0ND7S ,  C0ND7N , 
and  C0ND8 •  The  procedure  used  in  CONDI,  C0ND2 ,  C0ND3 ,  and 
C0ND5  is  that  of  Salacka  [Ref.  2:pp.  36-37]  modified  to 
account  for  contact  surfaces.  Appendix  C,  Fig.  C.5  is  a 
general  flowchart  for  these  four  routines.  Essentially,  an 
initial  estimate  is  made  of  the  characteristic  slopes,  ,  by 
enforcing  the  principle  of  domain  of  dependence.  The 
assumption  is  made  that  the  slopes  are  linear,  and  that  the 
characteristics  pass  through  point  B,  as  defined  in  Fig. 
3.7.  Then  q  and  A  are  computed  at  point  A,  which  in  turn 
yields  a  second  estimate  of  the  characteristic  slope,  ,\. 

The  two  slopes  for  each  characteristic  are  compared. 
If  they  are  not  equal  to  within  an  acceptable  error,  the  new 
estimate  of  \  is  used  to  repeat  the  process  until 
convergence  is  achieved.  All  three  characteristics  are 
handled  simultaneously.  The  value  of 

is  =  A-lt  (3.2) 

is  then  determined  for  each  characteristic  q+A,  q,  and  q-A. 


This  allows  linear  interpolation  of  Q,  R,  and  S  at  point  A. 
The  spatial  derivatives,  8<2/3s  and  9 y  as  are  computed  from 


t+dt  ◄ - ►  B 


Figure  3.7  Computational  Method  for  CONDI  Routine 

the  values  of  q  and  A  used  in  the  estimate  of  the  initial 
characteristic  slopes. 

The  "CONDI”  subroutine  is  used  when  there  are  no 
contact  surfaces  or  shocks  within  the  interval  on  either 
side  of  the  node,  and  flow  is  subsonic.  The  q+A  character¬ 
istic  uses  forward  differencing,  while  the  q-A  characteristic 
uses  a  backward  differencing  scheme  to  keep  the  initial 
domain  of  dependence  of  the  numerical  scheme  outside  the 
physical  domain  of  dependence. 

The  "C0ND2"  subroutine  is  a  backward  differencing 
algorithm  used  in  situations  when  one  or  more  discontinui¬ 
ties  are  in  the  right  interval,  or  if  flow  is  supersonic  to 
the  right.  Fig.  3.8  illustrates  the  computational  method 
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Figure  3.8  Computational  Method  for  C0ND2  Routine 

for  "C0ND2 " .  The  characteristic  slopes  of  q+A  and  q  are 
handled  as  in  "CONDI",  but  the  q-A  characteristic  slope  is 
determined  by  a  backward  rather  than  a  forward  scheme.  To 
interpolate  between  node  i  and  i+1  would  be  incorrect 
because  of  the  discontinuity  in  the  values  between  them. 
The  value  of  parameters  in  the  right  interval  are  interpo¬ 
lated  from  values  in  the  left  interval  by  assuming  the 
derivatives  do  not  change  between  the  intervals  up  to  the 
discontinuities  (i.e.,  a  shock  and/or  a  contact  surface). 

The  "C0ND3 "  subroutine  is  a  forward  differencing 
algorithm  used  in  situations  when  one  or  more  discontinui¬ 
ties  are  in  the  left  interval,  or  if  flow  is  supersonic 
traveling  to  the  left.  Figure  3.9  shows  the  computational 


discontinuity  B 


Figure  3.9  Computational  Method  for  C0ND3  Routine 

method  for  "C0ND3 " .  The  method  is  similar  to  that  in  C0ND2 
but  the  q-fA  characteristic  slope  has  to  be  interpolated  from 
values  of  node  i  instead  of  node  i-1  using  a  forward 
difference  scheme.  The  characteristic  slope  for  q-A  and  q 
are  estimated  by  forward  difference  schemes  using  values  of 
node  i+1  and  i  respectively. 

For  conditions,  such  as  that  illustrated  in  Fig. 
3.10,  where  a  discontinuity  is  in  the  interval  opposite  the 
direction  of  supersonic  flow  then  C0ND5  subroutine  is 
called.  A  mix  of  C0ND2  and  C0ND3  is  used.  For  flow  to  the 
right  the  method  of  C0ND2  is  implemented,  while  flow  to  the 
left  is  the  same  as  for  C0ND3 . 


discontinuity. 


q,  >  A  , 
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Figure  3.10  Example  Condition  for  "C0ND5"  Routine 

Subroutine  "C0ND4”  is  called  whenever  a  node  is  to 
be  crossed  by  a  discontinuity  from  either  direction.  The 
flowchart  for  "COND4"  is  shown  in  Appendix  C,  Fig.  C.6.  Two 
integer  vector  arrays,  LNODE  and  RNODE,  are  used  to  store 
the  following  information  on  the  node  being  crossed. 

A)  the  node  number,  I 

B)  the  value  of  SHOCK 

C)  the  value  of  CNTACT 

D)  the  current  value  of  the  integer  time  step,  J 

The  information  is  used  in  "CORRCT"  to  update  the  node  to 
the  next  time  level.  In  the  current  code  the  maximum  number 
of  nodes  that  can  be  crossed  during  any  time  step  is  two. 
Figure  3.11  shows  an  example  of  this  situation,  and  defines 
the  assignment  of  values.  Because  the  nodes  are  swept  from 
left  to  right,  the  left-most  node  is  assigned  to  LNODE, 
while  RNODE  applies  to  any  node  crossed  to  the  right  of  the 
LNODE  node.  The  values  of  ivT  are  set  to  zero,  so 
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Figure  3.11  Example  of  " COND4 "  Routine  Situation 


that  these  nodes  are  skipped  when  updating  occurs  in 
"SWEEP".  The  main  program  interprets  the  value  of  J  in 
LNODE  and  RNODE  to  determine  if  zero,  one,  or  two  nodes  need 
to  be  corrected  in  "CORRCT".  If  zero  nodes  are  crossed 
during  a  time  step  then  "CORRCT"  is  not  called. 

The  "C0ND6"  and  "C0ND8"  subroutines  are  designed  for 
future  extensions  of  the  entire  code.  Until  the  Rieman--: 
problem  illustrated  in  Fig.  3.6  is  coded,  then  any  situation 
after  the  shock  and  contact  surface  intersect  cannot  be 
computed.  However,  "SWEEP"  is  already  coded  to  call  "C0ND6" 
for  the  situation  illustrated  in  Fig.  3.12,  or  its  mirror 
image.  "C0ND8"  would  be  called  for  all  the  other  situations 
after  the  shock  and  contact  surface  interaction  has  taken 
place.  Currently,  both  subroutines  output  messages  on  the 
screen  describing  the  situation,  and  set  HALT  equal  to  1  to 
terminate  the  program. 


shock  contact  surface 


Figure  3.12  Example  of  "C0ND6"  Situation 

The  "C0ND7",  "C0ND7N",  and  "C0ND7S"  subroutines  are 
all  related  to  each  other.  They  are  intended  to  contain 
code  that  would  facilitate  the  shock  and  contact  surface 
intersection  and  solve  the  resulting  Riemann  problem.  The 
"C0ND7"  routine  is  called  when  the  contact  surface  and  shock 
moving  in  opposite  directions  would  cross  in  a  time  step 
within  an  interval.  The  routine  would  calculate  when  and 
where  within  the  time  and  space  domain  the  intersection 
would  occur  based  on  the  known  velocity  of  each  discontinu¬ 
ity  and  their  current  locations.  This  new  time  could  then 
be  used  to  rerun  "SWEEP"  with  the  shock  and  contact  surface 
exactly  on  top  of  each  other  at  the  end  of  the  new  time 
step.  The  "C0ND7N"  routine  is  for  the  -same  situation,  but 
when  a  node  is  crossed  by  either  discontinuity  during  the 
same  time  step.  The  same  procedure  is  done,  with  the 
requirement  to  ensure  the  node  crossed  is  properly  updated 


by  "CORRCT"  during  the  new  time  step.  The  "C0ND7S"  routine 
is  called  when  a  shock  and  contact  surface  are  located  at 
the  sane  point  at  a  time  other  than  zero.  The  subroutine 
would  contain  the  code  to  solve  the  Riemann  problem  set  up 
by  "C0ND7N"  or  "C0ND7 "  routines. 

6.  The  "CORRCT"  Routine 

Appendix  C,  Fig.  C.7  shows  the  flowchart  for 
"CORRCT".  The  "CORRCT"  routine  corrects  nodes  that  have 
been  crossed  by  a  discontinuity  (i.e.,  a  shock  or  a  contact 
surface)  or  are  straddled  by  both.  The  routine  takes  the 
information  from  LNODE  and  RNODE  arrays  to  determine  which 
node  or  nodes  are  to  be  updated,  and  if  there  is  any  shock 
and  contact  surface  interaction  at  the  nodes.  Then,  using 
subroutines  "DELTAX",  "SKJUMP",  "CSJUMP",  "EXTRAP", 
"INTERP",  and  "BBDRY"  the  concepts  from  Section  II  are 
imposed  to  calculate  the  jump  in  the  parameters  R,  Q,  S,  A 
and  q  at  the  node  in  question  to  the  new  time  level. 

7.  The  "BONDRY"  Routine 

The  "BONDRY"  routine  computes  the  new  values  of  the 
parameters  Q,  R,  and  S  at  the  boundary  nodes  for  time  step 
t.  The  concepts  of  Section  II  for  an  open  or  closed 
boundary  are  coded  as  shown  in  the  flowchart  of  "BONDRY" 
given  in  Appendix  C,  Fig.  C.8.  The  routine  uses  the 
"CONDI",  "C0ND2" ,  or  "C0ND3 "  routines  to  calculate  the 


correct  w,  and  then  solves  for  w  the  same  as  in  "SWEEP". 
If  a  shock  crosses  an  open  boundary  node  during  a  time  step. 


then  "SKJUMP"  routine  is  used  to  calculate  the  values  behind 


the  shock  at  the  node. 

8.  The  "SRFLCT"  Routine 

The  "SRFLCT"  subroutine  is  called  to  calculate  the 
values  of  Q,  R,  S,  q  and  A  at  the  boundary  node  when  a  shock 
reflects  from  a  solid  boundary.  Appendix  C,  Fig.  C.9  shows 
the  flowchart  of  "SRFLCT"  which  codes  the  analysis  in 
Section  II. D. 2  for  closed  boundary  shock  reflection.  The 
time  for  the  shock  to  reach  the  solid  boundary,  3tw,  is 
computed.  Then  the  excess  time  in  the  time  step,  3t,  is 
given  by 


3tex  *  3t  -  otw  (3.3) 

After  the  new  shock  velocity,  Vs, is  calculated,  then  the  new 
location  of  the  reflected  shock  is  known  from  multiplying  Vs 
by  ~J  fcex  and  adding  or  subtracting  from  the  boundary 
location. 

9.  The  "DELTAX"  Routine 

The  "DELTAX"  subroutine  is  used  in  "CORRCT"  and 
"BONDRY"  to  calculate  the  location  within  an  interval  of  a 
discontinuity  in  terms  of  x.  Figure  3.13  defines  the 
various  displacements  which  are  calculated. 

10.  The  11  SKJUMP"  Routine 


The  "SKJUMP"  subroutine  computes  the  conditions 
downstream  of  a  stationary  normal  shock  as  described  in 


Figure  3.13  Discontinuity  Location  within  an  Interval 


Section  II. D. 2.  The  "CORRCT",  "TRAK",  and  "BONDRY" 
subroutines  call  "SKJUMP"  when  required,  to  obtain  the  speed 
of  sound  (A),  entropy  (S)  ,  and  velocity  (q)  behind  the 
shock. 

11.  The  "CSJUMP"  Routine 

The  "CSJUMP"  subroutine  computes  the  velocity  and 
speed  of  sound  change  across  a  contract  surface  as  described 
in  Section  II. D. 3.  The  "CORRCT"  and  "TRAK"  routines  call 
"CSJUMP". 

12.  Numerical  Support  Routines 

The  "EXTRAP",  "INTERP",  and  "BBDRY"  routines  are 
used  by  "CORRCT",  "BONDRY",  "TRAK",  and  "DBURST"  to  extrapo¬ 
late  or  interpolate  data  to  the  face  of  a  discontinuity  or 
point.  In  particular,  "BBDRY"  is  an  interpolation  routine 
used  within  an  interval  near  a  boundary. 


The  four  output  subroutines  used  are  those  developed 
by  Salacka  [Ref.  2:pp.  39-40].  The  “BORDER"  and  "PLOT" 
routines  output  plots  of  pressure,  velocity,  density,  and 
modified  entropy  distributions  versus  a  non-dimens ion- 
alized  tube  length  at  preset  time  intervals.  "EXACT" 
creates  a  plot  of  the  exact  density  distribution  at  selected 
points  and  compares  it  with  the  computed  density 
distribution  for  a  selected  time  step.  The  "LIST"  routine 
outputs  to  files  8,  9,  and  10  tabular  listing  of  computed 
data.  The  files  are  sent  to  the  user's  permanent  storage 


disk. 


A  value  of  GRAPHS  equal  to  one  causes  "BORDER"  to 


be  called  in  the  main  program.  This  defines  plot  axis, 
labels,  and  headings.  "PLOT"  is  then  called  once  at  time 
zero,  and  at  every  time  step  set  by  SKIP. 

If  GRAPHS  is  set  equal  to  two,  then  "EXACT"  is 
called  every  SKIP  time  steps  to  plot  six  exact  density 
values  and  the  computed  density  distribution.  The  six  exact 
points  are: 

A)  The  two  boundary  points 

B)  The  head  and  tail  of  the  rarefaction  wave 

C)  .  A  point  just  behind  the  contact  surface 

D)  A  point  just  behind  the  shock. 


WWW 


The  location  of  the  exact  values  is  based  on  elapsed  time 
and  known  exact  values  for  wave  velocities  entered  with  the 
initial  conditions  for  the  problem. 

"LIST"  routine  is  called  every  SKIP  time  steps  when 
GRAPHS  is  equal  to  zero.  File  9  contains  a  tabular  listing 
of  the  Riemann  variables,  modified  entropy,  pressure, 
density  and  velocity  distributions,  elapsed  time,  discon¬ 
tinuity  velocity  and  location.  File  10  is  a  tabular  listing 
of  the  location  of  the  shock,  contact  surface,  and  the  head 
and  tail  of  the  expansion  wave  computed  at  every  SKIP  time 
steps.  File  8  contains  the  exact  values  of  the  location  of 
the  shock,  contact  surface,  and  the  head  and  tail  of  the 


expansion  wave. 
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Four  test  cases  were  run  using  the  E1DV2  code.  The 
purpose  was  to  verify  the  ability  of  the  code  to  determine 
the  unsteady  flow  process  and  correctly  simulate  varying 
boundary  conditions  and  flow  directions  in  the  Riemann  shock 
tube  problem.  The  shock  tube  is  illustrated  in  Fig.  4.1 
with  the  high  pressure  on  the  left.  The  tube  is  divided 
into  two  sections  by  a  diaphragm.  When  the  diaphragm  is 
burst,  the  pressure  equalizes  through  a  shock  wave  traveling 
into  the  expansion  chamber,  and  an  expansion  or  rarefaction 


diaphragm 


time  *  0 


P  Vexpanslon  wave  contact 

.  shock 


time  =  t 


Figure  4.1  Shock  Tube  at  t  =  0  and  t 


wave  traveling  into  the  compression  chamber.  A  contact 
discontinuity  is  behind  the  shock  wave,  traveling  at  the 
particle  velocity  [Ref.  6:p.  30], 

A.  TEST  CASE  1 

Test  Case  1  was  designed  to  demonstrate  shock  reflection 
and  expansion  wave  reflection  at  solid  boundaries.  Test 
Case  1  had  the  following  initial  and  boundary  conditions: 

1)  Pressure  and  density  ratios  equal  to  five,  with  high 
pressure  on  the  left 

2)  Temperature  ratio  equal  to  unity 

3)  Both  boundaries  were  closed 

4)  Diaphragm  was  located  at  x  =  0.5 

5)  Computational  mesh  had  101  nodes 

6)  Maximum  time  step,  JSTOP,  =  109,  with  SKIP  =  18. 

Plots  of  the  results  obtained  for  the  pressure,  density, 
velocity,  and  modified  entropy  distributions  are  shown  in 
Fig.  4.2  and  Fig.  4.3.  Fig.  4.2  shows  the  computation  up  to 
time  step  55,  while  Fig.  4.3  takes  the  computation  from  time 
step  56  to  termination.  The  output  of  "EXACT",  a  plot  of 
the  exact  density  compared  with  the  computed  density 
distribution  is  shown  in  Fig.  4.4.  Fig.  4.5  illustrates  the 
spatial  location  versus  time  for  exact  and  computed  shock, 
contact  surface,  and  head  and  tail  of  the  expansion  wave. 
Data  for  Fig.  4.5  were  taken  from  file  10  and  file  8.  Fig. 
4.6  shows  plots  of  pressure,  density,  modified  entropy,  and 
velocity  distributions  for  the  same  initial  and  boundary 
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Figure  4.2  Test  Case  1  (J  =  1  to  J  «  55) 
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Figure  4.6  Test  Case  1,  High  Pressure  on  Right 
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conditions  but  with  high  pressure  on  the  right.  Figure  4.6 
includes  data  corresponding  to  those  in  both  Fig.  4.2  and 
Fig.  4.3  for  the  high  pressure  initially  on  the  left. 


B.  TEST  CASE  2 

Test  Case  2  was  to  demonstrate  an  open  boundary, 
constant  pressure,  expansion  wave  interaction.  The  test 
case  was  run  with  the  following  initial  and  boundary 
condit ions : 


1,  Pressure  and  density  ratios  equal  to  5.3,  with  nigh 

pressure  on  the  left 

2)  Temperature  ratio  equal  to  unity 

’  The  left  ooundary  was  open,  with  a  constant  pressure 
ratio  across  the  boundary  of  4,5 

*  The  right  boundary  was  open,  with  an  adjustable  ores- 
sure  rat.o  t.nat  .n  effect  extended  the  .  e.ngt.n  or  the 
tube 

5  The  diaphragm  was  located  at  :<  =  ).~4 
•■>  Comput it iona  .  mesh  had  5  1  nodes 

Maximum  time  step .  23 TCP  *  : ~ .  JKI?  *  ?. 

Plots  of  the  results  for  the  pressure,  density,  velocity. 


mi  -ou.r.ed  m.trocy  list  r  .cut ,  :ns  ire  mown  .r,  ■  i  j.  . 

The  test  case  was  rerun  with  high  pressure  on  the  right. 
F.gure  4.3  snows  the  results  for  pressure,  density,  and 
temperature  -at  i  os  -.he  same.  The  i  i  iphragm  was  then  it  .<  = 
0.2f.  Boundary  conditions  were  reversed.  The  same  computa¬ 
tional  mesh,  and  output  control  were  used. 
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Figure  4.8  Test  Case  2,  High  Pressure  on  Right 
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C.  TEST  CASE  3 


Test  Case  3  was  designed  to  demonstrate  a  shock  exiting 
an  open  boundary  with  constant  pressure  maintained  at  the 
boundary  itself.  The  following  initial  and  boundary 
conditions  were  set: 

1)  Pressure  and  density  ratios  equal  to  5.5,  with  high 
pressure  on  the  left 

2)  Temperature  ratio  equal  to  unity 

3)  Left  boundary  was  closed 

4)  Right  boundary  was  open  with  a  constant  pressure  ratio 
across  the  boundary  of  unity 

5)  Diaphragm  was  located  at  x  ■  0.5 

6)  Computational  mesh  had  101  nodes 

7)  Maximum  time  step,  JSTOP,  »  73,  with  SKIP  *  18. 

Figure  4.9  shows  plots  of  the  results  obtained  for  pressure, 
density,  modified  entropy,  and  velocity  distributions. 
Similarly,  Fig.  4.10  shows  results  obtained  by  putting  the 
high  pressure  on  the  right,  reversing  the  boundary  condi¬ 
tions,  and  holding  everything  else  the  same. 

D.  TEST  CASE  4 

Test  Case  4  was  designed  to  demonstrate  a  lower  pressure 
ratio,  with  shock  wave  reflection.  The  initial  and  boundary 
conditions  were  set  as  follows: 

1)  Pressure  and  density  ratios  equal  to  3.2,  with  high 
pressure  on  the  left 

2)  Temperature  ratio  equal  to  unity 

3)  Both  boundaries  were  closed 
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4)  Diaphragm  was  located  at  x  *  0.5 

5)  Computational  mesh  had  901  nodes 

6)  Maximum  time  step,  JSTOP,  -  601,  with  SKIP  150. 

Plots  of  the  results  obtained  for  the  pressure,  density, 
modified  entropy,  and  velocity  distributions  are  given  in 
Fig.  4.11.  Figure  4.12  shows  the  results  obtained  with  the 
high  pressure  to  the  right,  SKIP  =  2  00,  and  holding  all 
other  parameters  constant. 


RCSSURL  VELOCITY 


SHOCK  TUBE  RESULTS 

FIRST  ORDER  N  -901 
DENSITY  RATIO  -  3.2  TEMP  RATIO 
PRESSURE  RATIO  -  3.2 


J. 00  0.25 


.00 

0.00  0.25 

J3  ~ 

>—• 

<  \  \ 
\  \ 


A.  RESULTS  OF  TEST  CASES 
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Overall,  in  the  four  specific  test  cases  which  were  run, 
the  expected  qualitative  flow  behavior  was  produced  by  the 
code.  Quantitatively,  either  experimental  data  or  further 
exact  solutions  are  needed  to  fully  validate  the  computa¬ 
tions.  The  results  of  each  test  case  will  be  discussed 
separately. 

Test  Case  1  results  in  Fig.  4.2  show  the  formation  of  a 
well-defined  shock  wave  traveling  to  the  right,  with  the 
contact  surface  crisply  defined  following  along  behind.  The 
entropy  drops  sharply  across  the  shock  and  remains  constant 
to  the  contact  surface  and  then  jumps  to  a  steady,  constant 
value  to  the  left  boundary  as  expected.  Pressure  and 
velocity  remain  constant  through  the  contact  surface. 
Velocity  is  positive  since  flow  is  to  the  right.  The 
expansion  wave  is  smeared  from  head  to  tail  over  the  correct 
range.  Figure  4.3  is  a  continuation  of  the  flow  problem 
with  a  sequential  display  of  the  results.  The  entropy 
continues  to  drop  as  the  reflected  shock  travels  back  toward 
the  contact  surface.  Notice  that  the  velocity  behind  the 
shock  is  zero.  At  the  left  boundary  as  the  head  of  the 
expansion  wave  reflects,  pressure  and  density  continues  to 
drop  while  velocity  is  zero  at  the  boundary.  The  density 
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distribution  shows  that  the  contact  surface  and  shock  are 
about  to  cross.  The  program  terminated  and  issued  a  message 
before  the  next  output  call  was  made,  when  the  code 
determined  that  the  contact  surface  and  shock  would  cross. 
The  exact  density  distribution  compared  with  the  ccrr  .-el 
density  distribution  in  Fig.  -i .  4  clear.'/  shows  the  she:-:  iru 
contact  surface  are  modeled  very  accurately  with  a  medium 
mesh.  The  expansion  process,  consisting  of  infinitesimal 
changes  propagating  along  the  character ist ics  is  fairly  well 
modeled.  Figure  4.5  shows  that  tracking  of  the  contact 
surface  and  shock  wave, in  comparison  with  exact  locations  at 
set  times,  is  excellent.  Using  q-*-  4  and  q--  based  on  condi¬ 
tions  behind  the  burst  diaphragm  as  estimates  of  the 
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[Ret.  9:pp.  181-191].  However,  an  examination  of  the 
tabulated  output  showed  that  the  solution  method  predicted 
changes  in  the  gradients  to  occur  at  similar  locations  to 
those  computed  using  the  q-*-4  and  q-4  estimates.  Consequent¬ 
ly,  while  there  may  be  some  inaccuracy  in  the  numerical 
solution,  the  method  of  tracking  the  head  and  tail  of  the 
expansion  wave  is  acceptable. 

The  ability  to  reproduce  the  same  conditions  but  with 
flow  to  the  left  is  demonstrated  clearly  in  Fig.  4.6.  The 


•ntropy,  density,  velocity,  and  pressure  distributions  are 
mirror- images  of  those  for  the  case  of  high  pressure  on  the 
left. 

Test  Case  2  results  in  Fig.  4.7  show  that  at  time  zero 
with  a  pressure  ratio  at  the  boundary  of  4,  which  is  lower 
than  the  preset  value  inside  the  tube  of  5,  an  expansion 
wave  traveling  inwards  is  generated.  Outflow  is  seen  in  the 
negative  velocity  distribution  developing  near  the  left 
boundary,  where  negative  velocity  implies  flow  traveling  to 
the  left.  The  contact  surface  and  shock  are  clearly  formed 
to  the  right  of  the  off-center  diaphragm,  and  travel  to  the 
right.  The  expansion  wave  generated  at  the  diaphragm 
travels  to  the  left  and  intersects  with  the  expansion  wave 
generated  at  the  left  boundary.  The  pressure  distribution 
shows  the  propagation  of  these  waves  as  the  pressure  drops 
behind  the  head  of  each  wave.  With  conditions  reversed,  so 
that  the  high  pressure  is  to  the  right  of  the  diaphragm,  tne 
results  in  Fig.  4.3  are  a  mirror-image  with  velocity 
negative  for  flow  to  the  left,  and  positive  for  the  outflow 
to  the  right. 

Test  Case  3  results  in  Fig.  4.9  demonstrate  tha*  * 
shock  wave  meeting  an  open  boundary  where  the  press'.-'' 
held  constant,  is  correctly  computed  to  exit  the  tu.-p 
density  distribution  shows  the  jump  increase  a  - 
shock,  then  another  jump  increase  across  ft* 
surface.  Then  density  plunges  down  as  \r  .  . 
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travels  back  into  the  tube  when  the  shock  exits.  The 
pressure  ratio  at  the  open  right  boundary  was  held  at  unity, 
as  seen  in  the  pressure  distribution  which,  at  the  boundary 
behaves  similarly  to  the  density  distribution.  The  pressure 
at  the  boundary  would  be  required  to  increase  as  the  shock 
passed  by,  however  the  enforcement  of  the  constant  pressure 
locally  causes  an  expansion  wave  to  form  traveling  to  the 
left.  The  entropy  remains  constant,  at  the  value  behind  the 
shock,  from  the  boundary  back  to  the  contact  surface.  There 
is  a  jump  in  entropy  across  the  contact  surface.  The 
velocity  distribution  shows  that  as  the  shock  travels  to  the 
right,  the  velocity  jumps  up.  When  the  expansion  wave 
generated  at  the  right  boundary  travels  inward  the  velocity 
continues  to  increase  in  magnitude  while  flowing  to  the 
right.  The  expansion  wave  generated  at  the  diaphragm  can  be 
seen  traveling  to  the  left  in  the  plots  of  pressure, 
density,  and  velocity.  Reversing  the  situation  and 
computing  conditions  for  flow  to  the  left,  in  Fig.  4.10, 
results  in  a  mirror-image  in  the  pressure,  density,  and 
entropy  distributions.  The  velocity  becomes  negative  since 
flow  is  to  the  left. 

Test  Case  4  is  a  repeat  of  the  first  test  case  but  with 
a  pressure  ratio  of  3.2  and  a  fine  computational  mesh.  In 
Fig.  4.11  the  density  steps  up  across  the  shock  and  contact 
surface  as  they  travel  to  the  right.  When  the  shock 
reflects  the  density  jumps  again.  The  velocity  distribution 


shows  that  the  velocity  drops  to  zero  behind  the  reflected 
shock,  after  it  had  jumped  up  across  the  shock  when  it  was 
traveling  to  the  right.  The  instability  seen  in  the 
velocity  and  pressure  distributions  near  the  midpoint 
occurred  during  the  first  100  time  steps.  A  numerical 
instability  appears  to  generate  at  the  diaphragm  at  time 
zero  for  pressure  ratios  less  than  about  3.0  which  can  be 
severe.  For  the  conditions  in  this  test  case  the  transient 
instability  damped  out  after  the  first  100  time  steps.  The 
shock,  contact  surface,  and  expansion  wave  are  nevertheless 
seen  to  be  accurately  modeled.  Figure  4.12  demonstrates 
that  reversed  conditions  result  in  mirror-images  of  the 
computed  conditions. 

B.  CURRENT  LIMITATIONS  OF  THE  CODE  E1DV2 

The  current  limitations  of  the  E1DV2  code  for  modeling 
quasi-one-dimensional  inviscid  flow  are  in  two  categories. 
First,  there  are  numerical  limitations  in  obtaining 
solutions  with  the  present  code  for  different  initial  and 
boundary  conditions.  Second,  the  present  coding  limits  what 
flow  situations  can  be  treated. 

The  numerical  instability  which  occurs  at  low  pressure 
ratios  during  the  first  few  time  steps  can  be  severe.  The 
assumption  in  the  current  code  is  that  the  shock  forms 
immediately,  which  at  high  pressure  ratios  does  not  pose  any 
problems.  The  mathematical  modeling  of  the  bursting 
diaphragm  in  subroutine  "D BURST"  adequately  reduces  the 


transient  instability  for  high  pressure  ratios,  but  not  for 
low  pressure  ratios. 

A  second  rather  different  numerical  limitation  was 
identified  in  numerically  detecting  the  location  of  the  head 
and  tail  of  the  expansion  process.  Because  the  expansion 
process  is  not  a  sharply  defined  front,  fixing  the  precise 
locations  of  the  head  and  tail  waves  numerically  at  a  given 
time  is  difficult.  However,  the  method  currently  used  does 
accurately  track  the  computed  propagation  of  the  wave  along 
the  characteristics.  The  characteristics  are  modeled 
currently  as  straight  lines.  A  higher  order  curve  would 
describe  them  more  accurately. 

The  current  code  is  limited  to  tracking  two  discontinui¬ 
ties  and  the  expansion  wave  (i.e.,  a  shock  and  a  contact 
surface) .  Thus  any  situation  which  would  generate  another 
discontinuity  can  not  be  computed.  The  code  can  determine 
when  this  will  occur,  and  will  issue  a  message  explaining 
the  situation  before  terminating.  Thus  the  shock  colliding 
with  a  contact  surface  will  currently  terminate  the  program. 
A  boundary  pressure  that  is  higher  than  the  pressure  inside 
the  tube  would  also  create  a  compression  wave  or  possibly  a 
shock  wave.  This  condition  will  also  terminate  the  program. 
Simulation  of  the  process  of  a  shock  forming  over  an 
extended  distance  and  time  as  compression  waves  pile  on  top 
of  each  other  and  strengthen,  would  also  require  additional 


code. 


Finally,  the  variation  of  gamma  that  occurs  across  the 
contact  surface  in  a  wave  rotor  cannot  be  handled  currently 
since  E1DV2  assumes  a  single  constant  gamma  throughout. 


Towards  the  development  of  a  one-dimensional  code  for 
wave  rotor  applications,  methods  for  tracking  and  correcting 
conditions  across  a  contact  discontinuity,  applying  open  and 
closed-end  boundary  conditions,  accounting  for  shock  wave 
and  contact  surface  interaction  were  devised  and  were 
presented  here  in  detail.  The  EULER1  Fortran  Code  [Ref.  2] 
was  revised  to  become  the  E1DV2  code  with  the  following 
additional  capabilities: 

1)  tracking  of  the  contact  surface  and  expansion  wave 

2)  imposing  high  pressure  initially  on  the  right  side  of 
the  diaphragm 

3)  correct  jump  conditions  across  the  contact  surface 

4)  allow  open  boundary  conditions  with  constant  pressure 
specified  and  allow  exiting  of  shock  and  expansion 
waves 

5)  allow  closed  boundary  conditions  and  model  shock  and 
expansion  wave  reflections 

6)  improved  numerical  accuracy  from  time  zero  to  the 
maximum  time  step. 

The  code  was  tested  on  the  shock  tube  problem  under  four 
different  initial  and  boundary  conditions  with  excellent 
results . 

The  following  extensions  are  recommended  to  make  the 
code  suitable  for  wave  rotor  applications: 

1)  Solve  the  Riemann  problem  at  the  momerft  when  the  shock 
and  contact  surface  intersect 


2)  Add  additional  code  to  track  more  than  two  discontinu¬ 
ities  and  one  expansion  wave 

3)  Incorporate  a  variable  value  of  gamma  into  the  code 

4)  Update  the  characteristic  curve  approximation  from 
linear  to  a  higher  order  polynomial 

5)  Add  code  necessary  to  handle  open  boundary  conditions 
with  inflow 

6)  Improve  on  the  numerical  computation  at  time  zero  for 
low  pressure  ratios  across  the  diaphragm 

7)  Enable  boundary  conditions  to  be  variable  during 
program  execution  to  allow  wave  rotor  cycle  design 

8)  Compare  computations  using  the  code  with  experimental 
data  where  available 

9)  Extend  the  code  to  two  dimensions. 

These  extensions  may  require  various  degrees  of  effort. 
However,  the  ability  of  the  QAZ1D  solution  method  to  be 
extended  rather  simply  to  describe  viscous  multi-dimensional 
flow  suggests  that  such  efforts  would  be  justified. 


APPENDIX  A 


DERIVATIONS  OF  EQUATIONS 

LIST  OF  VARIABLES 

Speed  of  sound 

Specific  heat  at  constant  pressure 
Specific  heat  at  constant  volume 
specific  internal  energy 
Specific  enthalpy 
Static  pressure 
Modified  Riemann  variable 
Reversible  heat  transferred 
Velocity  magnitude 
Modified  Riemann  variable 
Gas  constant 

Modified  entropy  (non-dimensional  where  S  =  S 

Static  temperature 

Time 

Velocity  relative  to  a  standing  shock  wave 
Specific  volume 

Incoming  Mach  Number  relative  to  a  stationary 
shock  wave 

Work 

Density 

Ratio  of  specific  heats 


B.  DERIVATION  OF  SHOCK  JUMP  EQUATION  FOR  HIGH  PRESSURE  ON 

THE  RIGHT 

The  analytical  equation  relating  the  non-dimensionalized 
extended  Riemann  variable,  R,  change  through  a  normal  shock 
is  derived  below  similar  to  that  for  the  Q  variable  change 
when  high  pressure  is  on  the  left  [Ref.  Z  :Appendix  A]. 

Figure  A.l  shows  a  shock  moving  to  the  left  with 
velocity  Vs.  Subscript  A  is  always  associated  with 
parameters  to  the  left  of  any  discontinuity,  and  B  with 
those  on  the  right.  To  enable  normal  shock  relations  to  be 
used  the  system  requires  a  velocity  in  the  opposite  direc¬ 
tion  but  of  equal  magnitude  be  imposed  upon  it.  The 
coordinate  system  was  defined  such  that  vectors,  such  as 
velocity,  directed  to  the  right  are  positive  while  those  to 
the  left  are  negative. 
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Figure  A.l  Shock  Wave  with  High  Pressure  on  Right 


Since  relative  incoming  Mach  Number,  w,  relates  the 
velocity  downstream,  or  that  region  into  which  the  shock  is 
traveling,  to  the  speed  of  sound  in  that  region  [Ref.  10:pp. 


114-154]  then 


w 


q  -V 
^A  s 


(Al) 


w  is  a  positive  value  because  qA,  though  negative,  is  small¬ 
er  in  magnitude  than  the  equally  negative  shock  velocity, 
Vs.  The  speed  of  sound.  A,  is  positive  by  definition. 

The  appropriate  extended  Riemann  variable  in  this  case 
is  R,  defined  as 

R  =  q  -  AS  (A2 ) 


Since  q  is  negative,  this  can  be  rewritten  as 


R  =  ~  ( I  ql  +  AS) 


(A3) 


f 

to  emphasize  that  R  is  the  Riemann  variable  associated  with 
the  lesser  change  across  the  shock  [Ref.  3:pp.  18-21].  We 
adopt  the  pattern  of  non-dimensionalizing  velocity  by  the 
speed  of  sound,  pressure  and  density  by  corresponding  values 
downstream  of  the  shock.  Using  the  entropy  S,  defined  by 
Verhoff  [Ref.  l:p.  2]  as 


then 
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The  ratios  of  pressure,  density,  and  sonic  velocity  across  a 
normal  shock  from  Zucker  [Ref.  7:p.  151]  and  Shapiro  [Ref. 
10 :p.  118]  in  terms  of  Mach  Number  are 
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The  first  term  on  the  right  side  of  equation  (A5)  can  be  ex¬ 
pressed  in  terms  of  w  by  applying  simple  continuity  in  mass. 
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Taking  the  areas  as  equa],  this  becomes 
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Combining  equations  (A6) ,  (A7) ,  (A8)  ,  and  ( A 1 3 ; 

equation  (A  5)  gives 
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C.  DERIVATION  OF  CONTACT  SURFACE  JUMP  EQUATION 

The  analytical  egression  for  the  ratio  of 
velocity  across  a  constant  surface  is  derived  -  n 
The  first  law  of  thermodynamics  is  stated  as 


work  done 
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surround  i  r.a 
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change  of 
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+ 

or,  in  differential  form 


de  —  9 qp  +  9W 

Internal  energy  is  defined  as 


(A15; 


e  =  h  -  pv 


(A16) 


then 
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de 


dh  -  pdv  -  vdp 


(A17) 


The  incremental  work  is  given  by 

3w  =  pdv  ( A18 ) 

The  heat  addition  is  determined  from  the  definition  of  modi¬ 
fied  entropy. 
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Substituting  equations 
(A15)  yields 
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(A17) ,  (A18) ,  and  (A19) 


(A19) 

into  equation 


dh  -  pdv  -  vdp  =  -yTdS  -  pdv 


Canceling  like  terms,  and  rearranging  gives 

dp  =  X  TdS  +  (— )dh  (A20) 

v  v 

Enthalpy,  h,  is  defined  as 


h  =  CpT 


(A21) 


For  a  perfect  gas  the  sonic  velocity  is  given  by 

A  -  (yRgT)V2 

thus 

dA  -  jfr  RgT)  1/2  (dT/T)  (A22) 

Substituting  equation  (A21)  and  then  (A22)  into  equation 
(A20)  ,  equation  (A.20)  becomes 

dp  -  ±TdS  +  (— )  d  ( CpT) 

-  ^TdS  +  (^)CpdT 

v  —  1  2dA 

-  *TdS  +  (“)  CpT  (~^~)  (A23) 

For  an  ideal  gas, 

Y  *  Cp/  Cy 

and 

Cp  =  Cy  +  R  Q 


Combining  these  two  equations,  and  rearranging 
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Substituting  equation  (A24)  into  equation  (A23) ,  and 
observing  that  pressure  remains  constant  through  the  contact 
surface 
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Eliminating  T/v,  equation  (A25)  becomes 
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Then  becomes  a  non-dimensionalized  entropy.  Using  the 
notation  of  S  =»  S/Rq  now  for  the  non-dimensionalized  form, 
so  that  by  integrating  both  sides 
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which  can  be  written  as 


Aa/Ab  =  exp 
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Definitions  used  in  this  development,  excluding  modified 
entropy  were  from  [Ref.  ll:pp.  87-125J . 


D.  DERIVATION  OF  EXTENDED  RIEMANN  VARIABLE  CHANGE  ACROSS  A 
CONTACT  SURFACE 

An  analytical  expression  for  the  change  in  the  non- 
dimensionalized  extended  Riemann  variable,  Q,  across  a 
contact  surface  traveling  right  is  derived  here.  Figure  A. 2 
depicts  a  contact  surface  moving  right  with  velocity  q. 


Entropy 


Figure  A. 2  Contact  Surface  Traveling  Right 
Subscript  Notation 


Values  of  parameters  to  the  left  of  the  interface  are 
denoted  by  subscript  A',  and  those  to  the  right  with  sub¬ 
script  B'.  As  illustrated,  parameters  to  the  left  of  the 
shock  have  subscript  A,  and  to  the  right  subscript  B.  All 
velocities  are  non-dimensionalized  by  sonic  velocity 
immediately  downstream  of  the  discontinuity.  Thus  using  the 
definition  of  the  extended  Riemann  variable, 


I] 


q  +  AS 


(A28) 


then 


°A'  ”9b*  qA*  +AA,SA'  qB'  +AB,SB' 


1 


^A'  -<^B'  AA' 

*SB' 


(A29) 


Since  velocity  is  constant  through  a  contact  surface,  so 
that 
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then  equation  (A29)  becomes 
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Using  equation  (A27) ,  this  can  be  written 
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Multiply  each  side  by  ^g./Ag,  and  since  AR|  =  Aa  then 
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For  a  contact  surface  traveling  to  the  left,  as 
illustrated  in  Fig.  A. 3,  the  derivation  is  entirely  similar  but 
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Figure  A. 3  Contact  Surface  Traveling  Left 
Subscript  Notation 


using  the  R,  extended  Rieinann  variable  with  the  result 
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E.  DERIVATION  OF  OPEN  BOUNDARY  CONDITIONS  WITH  CONSTANT 
PRESSURE 

The  equation  for  density  and  temperature  ratios  in  terms 
of  pressure  and  entropy  are  derived  below.  The  definition 
of  modified  entropy  in  non-dimensional  form  is 


s  ■  FT  -  ln 


(A33) 
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where  P  and  p  are  non-dimensional ized  ratios. 


(S 


“  ^zr>  (“Y  (Y“l ) )  *  In  (P/pY ) 


or 


f  <S  -^j)  (-y(Y-D)J 


exp  1  ~  -  p/pY 


Rearranging, 


thus 


y  ,  ,  [(s‘^r)("Y(Y"1,)J 

PT  «  [Ptexp  Y  1  J] 


1A  (S'^>(“Y(Y“1>>  -1A 
P  «  Pi/Y(exp  Y  1  ) 


p  *  {^{exp 


1  (S y-)  (-y(y-l) )  -1/y 


)} 


Using  the  perfect  gas  law* 


fl#  «-4)H(ri))  -1/y 

0  '  ^<«*P  Y  )) 


-Y  1  r*)  (-Y(y-l) ) 

Y  -  —(exp  Y  1 


pT 


multiply  each  side  by  T  and  divide  by  o"Y  then 


y  (S (-y(y-l)) 
T  *  o yo  (exp  y  ) 


thus 


T  =  pT  [exp  y  i  ] 


(A35) 
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APPENDIX  B 


OPERATION  OF  E1DV2  ON  THE  NPS  VM/CMS  SYSTEM 

A.  TERMINAL  REQUIREMENTS 

Terminal  requirements  depend  on  the  desired  output.  Any 
terminal,  such  as  the  IBM  3278,  connected  into  the  VM/CMS 
system  is  adequate  for  a  tabular  listing  or  Disspla  Metafile 
of  data.  The  Disspla  Metafile  stores  graphical  data  for 
display  using  DISSPOP  commands.  If  output  of  the  plots  on  a 
monitor  screen  is  desired,  an  IBM  3277-TEK618  dual  screen 
terminal  must  be  used.  Table  B.I  lists  terminal  require¬ 
ments  for  different  outputs. 


TABLE  B.I 

TERMINAL  AND  OUTPUT 


OUTPUT 


TERMINAL 


Tabular  listing  of  data  Any  terminal 

Graphical  data  in  format 

for  VERSATEC  printing  Any  terminal 

Graphical  data  plotted 

on  screen  IBM  3277-TEK618 

B.  PROCEDURE 

1.  Editing  Variables  and  Program  Setup 

To  change  the  initial  and  boundary  conditions, 
graphical  output,  and  grid  size  to  run  the  program  under 


different  conditions,  the  following  must  be  done.  First, 
pick  the  appropriate  terminal  from  Table  B.I  and  log  on. 
Enter  the  editor  by  issuing  the  command 

x  E1DV2  FORTRAN  A 

Table  B. II  lists  the  code  variables  that  will  require 
editing  and  their  meaning.  In  addition,  to  these  variables 
if  graphical  plots  are  to  be  outputed  then  enter  the 
"BORDER”  and  "EXACT"  subroutines  and  edit  the  lines: 

"FIRST  ORDER  N  =  ???" 

"DENSITY  RATIO  =  ???  TEMP  RATIO  =  ???» 

"PRESSURE  RATIO  =  ???" 

To  issue  the  output  graphs  to  the  screen  comment  out  the 
lines  "COMPRS",  and  use 

CALL  TEK618 

Otherwise,  if  a  Disspla  Metafile  of  the  graphics  plot  is 
desired,  comment  out  the  "TEK618"  line,  and  use 

CALL  COMPRS 

These  lines  are  in  the  "Set  up  graphics  plot  of  variable" 
loop  in  the  main  program. 

To  store  these  changes,  hit  the  enter  key,  and  type 
"file".  Select  the  enter  key  once  more.  The  program  is  now 
ready  to  be  compiled  and  executed. 

2 .  Commands 

To  compile  and  execute  the  E1DV2  code  on  the  VM/CMS 
system  perform  the  following  commands  exactly  as  typed: 


TABLE  B. II 


Variable 

DIMENSION  statement 
for  arrays  A(N) , , 
XARRAY  (N) 

N 

GRAPHS 

SKIP 

JSTOP 

TRI 

PRI 

DRI 

QLI 

QRI 

LBDPRI 

LBDTRI 


EDITABLE  VARIABLES 
Edit 

set  the  dimension  of  the  arrays 
equal  to  the  number  of  nodes 
in  the  grid 

set  to  the  number  of  nodes  in  the 
grid 

set  to:  0  for  tabular  listing  of 
data 

1  for  plot  of  density, 
entropy,  pressure  and 
velocity  distributions 

2  for  plot  of  exact  density 
distribution  compared  to 
computed  density 
distribution 

set  to  number  of  time  steps  between 
calls  to  output  routines 

set  to  maximum  number  of  time  steps 

set  to  initial  temperature  ratio 
across  diaphragm 

set  to  initial  pressure  ratio  across 
diaphragm 

set  to  initial  density  ratio  across 
diaphragm 

set  to  initial  velocity  left  of 
diaphragm 

set  to  initial  velocity  right  of 
diaphragm 

set  to  initial  pressure  ratio 
across  left  boundary 

set  to  initial  temperature  ratio 
across  left  boundary 


LBDORI 


set  to  initial  density  ratio  across 
left  boundary 


TABLE  B. II  (CONTINUED) 


RBDDRI 

RBDPRI 

RBDTRI 

G 


set  to  initial  density  ratio  across 
right  boundary 

set  to  initial  pressure  ratio  across 
right  boundary 

set  to  initial  temperature  ratio 
across  right  boundary 

set  to  desired  value  of  gamma 


EE 


LWPRES 


LBNDRY 

RBNDRY 

LBDPRS 

RBDPRS 


XINIT 


set  to  desired  error  tolerance  for 
calculation  of  characteristic  slope 
(i.e.,  0.1D-8) 

set  to:  2  for  low  pressure  on  right 
side  of  diaphragm 
3  for  low  pressure  on  left 
side 

set  to:  1  for  closed  boundary 
0  for  open  boundary 

set  to:  0  for  constant  pressure  at 
left  boundary 

1  for  pressure  that  adjusts 
at  the  left  boundary 

set  to  location  of  diaphragm 


VHEAD 

VTAIL 

VCDE 


set  to  exact  velocity 
expansion  wave 

for 

head 

of 

set  to  exact  velocity 
expansion  wave 

for 

tail 

of 

set  to  exact  velocity 

for 

contact 

surface 

V3E 


set  to  exact  velocity  for  shock 
wave 


DLCD 


set  to  exact  density  behind  contact 
surface 


DLSH 


set  to  exact  density  behind  shock 


TABLE  B. II  (CONTINUED) 


SIGMA (1,2)  set  to  diaphragm  location  (i.e., 

SIGMA (2 , 2 )  0.5D00) 

SIGMA (3, 2) 

SIGMA (4 ,2) 

Y  There  are  four  cases  where  Y 

appears  in  the  program 
edit  as  follows: 

Comment  out  Y  =  (Integer#) ,  use: 
first,  Y  =  (N+l)/2  if  LWPRES  =  2 
second,  Y  =  (N+3)/2 
and 

third,  Y  =  (N-l)/2  if  LWPRES  =  3 
fourth,  Y  =  (N+l)/2 
for  diaphragm  at  0.5 
Comment  out  those  above,  if  dia¬ 
phragm  at  another  node.  Set 
first,  Y  =  (Integer  #)  of  node  for 
LWPRES  =  2 
second,  Y  =  (#  +  l) 
and 

third,  Y  =  (Integer  #  -  1)  of  node 
for  LWPRES  =  3 
fourth,  Y  =  (#) 


1)  Increase  the  virtual  memory  by  entering 

DEFINE  STORAGE  1M 

2)  To  return  to  CMS  environment  enter 

I  CMS 

3)  To  compile  the  program  enter 

FORTVS  E1DV2 

The  screen  will  display  messages  as  it  compiles  each 
routine  and  when  finished  a  ready  symbol  appears. 

4)  To  execute  the  program,  enter 


DISSPLA  E1DV2 


The  message 


...YOUR  FORTRAN  PROGRAM  IS  NOW  BEING  LOADED... 

...EXECUTION  WILL  SOON  FOLLOW... 

should  appear,  followed  by 

. . . EXECUTION  BEGINS . . . 

If  at  a  TEK618  terminal  with  GRAPHS  equal  to  1  or  2 
then  the  screen  on  the  TEK618  will  begin  plotting  the 
selected  graph.  A 

...press  ENTER  to  continue... 

message  will  appear  on  the  3277  terminal.  If  a  copy  of  the 
plot  is  desired,  do  so  now  before  pressing  the  enter  key. 
After  pressing  the  enter  key  on  the  3277  terminal,  the  plot 
will  be  erased  and  the  program  will  terminate.  Proper 
termination  will  result  in 

END  OF  DISSPLAY  9.2  ####  VECTORS  IN  1  PLOT... 
appearing,  followed  by  a  ready  message. 


If  GRAPHS  was  set  to  0,  then  proper  termination 
would  be  a  ready  message.  The  tabular  listing  of  pressure, 
density,  velocity,  entropy,  and  Riemann  variables  will  be  in 
"FILE  FT09F001."  The  exact  location  of  the  shock,  contact 
surface,  and  expansion  wave  with  elapsed  time  will  be  in 
"FILE  FT08F001."  The  computed  location  of  the  shock, 
contact  surface,  and  expansion  wave  with  elapsed  time  will 
be  in  "FILE  FT10F001." 

A  second  method  to  compile  and  execute  the  program, 
plus  provide  the  files  with  a  name  is  to  create  the  follow¬ 
ing  EXEC  file  on  the  user's  disk. 

FI  9  DISK  FILE09  LISTING  A (PERM 
FI  10  DISK  FILE10  LISTING  A (PERM 
FI  8  DISK  FILE08  LISTING  A (PERM 
FORTVS  &1 

A  possible  file  name  and  required  file  type  would  be 

RUN  EXEC 

To  compile  the  program  and  define  the  output  files,  enter 

.  RUN  E1DV2 

After  compiling  is  finished,  and  the  ready  message  appears, 
enter 

DISSPLA  E1DV2 

The  program  executes  as  outlined  before. 


APPENDIX  C 


FLOWCHARTS 

Flowcharts  for  major  routines  in  E1DV2  are  shown. 
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set  initial  and 


boundary  conditions 


Figure  C.l  Main  Program  Flowchart 


Figure  C.2  "DBURST"  SubrovJtine  Flowchart  . 
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Figure  C.2  (Continued) 
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Figure  C.3  (Continued) 
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(for  left  boundary  node) 

Figure  C.4  "SWEEP”  Subroutine  Flowchart 
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Figure  C.4  (Continued) 


Figure  C.4  (Continued) 
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Figure  C.5  General  "CONDI,  2,  3,  and  5"  Subroutine  Flowchart 
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Figure  C.6  "CONDV1  Subroutine  Flowchart 
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Figure  C.7  (Continued) 
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Figure  C.8  "BONDRY"  Subroutine  Flowchart 
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APPENDIX  D 


(UlM-10 

VERSION  2 
( E10V2 ) 

THIS  PROGRAM  SOLVES  THE  EULER  EQUATIONS 
EXPRESSEO  IN  A  QUA2I-0NE  DIMENSIONAL 
STREAMLINE  COORDINATE  SYSTEM. 

AUTHOR  -  LT.  O.T.  JOHNSTON, FEB  1987 


BASED  ON  THE  EULER1  COOE  BY 
T.F.  SALACKA ,  DEC  1985 

FEATURES  OF  THIS  VERSION  12) 

♦  ORDER  OF  SPATIAL  DERIVITIVES 

♦  NUMBER  OF  SPATIAL  DIMENSIONS 

♦  DISCONTINUITIES  TREATED! 

SHOCKS 

CONTACT  DISCONTINUITIES 
EXPANSION  HAVES 

♦  HIGH  PRESSURE  SIDE 

LEFT 

RIGHT 


FIRST 

ONE 


►♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

►  ♦ 

>  CONVENTIONS  ANO  DEFINITIONS  ♦ 

»  ♦ 

*♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦•♦♦♦♦♦♦♦ 

— NON- DIMENSIONING  CONVENTION - 

ALL  VELOCITIES  NON-DIM.  BY  THE  SOUND  SPEED  ON 
THE  LOH  PRESSURE  SIDE  OF  THE  DIAPHRAGM. 

ALL  PRESSURES,  DENSITIES,  ANO  TEMPERATURES  ARE 
NON-OIM.  8Y  THEIR  INITIAL  VALUES  ON  THE  LOW 
PRESSURE  SIDE  OF  THE  DIAPHRAGM. 

SPACIAL  OISTANCE  IS  NON-OIM. BY  OVERALL  LENGTH. 
ENTROPY  IS  NON-OIM.  BY  THE  GAS  CONSTANT,  R. 

TIME  IS  NON-OIM.  BY  UENGTH/SOUNO  SPEED). 
VELOCITIES  ANO  DISTANCES  ARE  OEFINED  POSITIVE 
TO  THE  RIGHT. 

.........  SUBSCRIPT  NOTATION— - - - 


I  -  SPACIAL  NODE  « 1  TO  N  FROM  LEFT  TO  RIGHT) 

J  -  TIME  LEVEL  )0  IS  THE  INITIAL  CONDITION) 

X  -  DENOTES  WHICH  CHARACTERISTIC  HAVE  IS  BEING 
DEALT-  HITHi 

1  »  Q*A  2  •  Q  5  ■  G-A 

L  -  DENOTES  WHICH  TYPE  OF  DISCONTINUITY  IS 
BEING  DEALT  WITH: 

1  ■  SHOCK  2  ■  CONTACT  DISCONTINUITY 
5  ■  HEAO  OF  EXPANSION  HAVE 
A  ■  TAIL  OF  EXPANSION  HAVE 


c 

C  A  -  speed  of  sound 

C  *A  -  DENOTES  THE  VALUE  OF  A  VARIABLE  AT  THE  NODE 

C  TO  THE  LEFT  OF  A  DISCONTINUITY .  *  CAN 

C  BE  ANY  VARIABLE  NAME. 

C  AR  -  THE  RATIO  OF  SOUNO  SPEED  ACROSS  A 

C  SHOCK,  A/BI  SHOCK  MOVING  RIGHT  )  ,B/AI SHOCK  MOVING  LEFT) 

C  «B  DENOTES  THE  VALUE  OF  A  VARIABLE  AT  THE  NODE 

C  TO  THE  RIGHT  OF  A  DISCONTINUITY. 

C  BORY  -  3  DENOTES  LEFT  BOUNDARY, 2  THE  RIGHT  BOUNOARY 

C  COUNT  -  COUNTER  FOR  GRAPHICS  ROUTINES 

C  OARRAY  -  ARRAY  OF  DENSITY  FOR  PLOTTING 

C  OELT  -  TIME  STEP 

C  OLCO  -  OENSITY  BEHIND  THE  CONTACT 

C  DISCONTINUITY  IN  THE  EXACT  SOLUTION. 

C  OLSH  -  OENSITY  BEHIND  THE  SHOCK  IN  THE 

C  EXACT  SOLUTION. 

C  DQ  -  THE  JUMP  IN  VELOCITY  ACROSS  THE  SHOCK 
C  DIVIDED  BY  THE  SOUNO  SPEED  AT  B) RIGHT)  OR  A( LEFT ) 

C  ORI  -  INITIAL  DENSITY  RATIO  ACROSS  THE  SHOCK 

C  EE  -  DESIRED  PRECISION  FOR  CHARACTERISTIC  CALCULATIONS 

C  G  GAMMA  (RATIO  OF  SPECIFIC  HEATS) 

C  GRAPHS  -  FOR  GRAPHICAL  OUTPUT,  0»N0NE  (TABULAR) 

C  1 “PLOTS  ALL  VARIABLES 

C  2 “COMP ARES  DENSITY  NITH  EXACT  SOLUTION 

C  G1  -  l/IG-l) 

C  G2  2/IG-l) 

C  H  l/( N-l ) 

C  HALT  -  TERMINATES  PROGRAM  IF  l.SET  BY  CONDITIONS  NOT  CODED 

C  12  NUMBER  OF  THE  NOOE  TO  THE  RIGHT  OF  A 

C  DISCONTINUITY. 

C  JSTOP  -  NUMBER  OF  TIME  LEVELS  TO  BE  CALCULATED 

C  LBOOR  -  LEFT  BOUNDARY  OENSITY  RATIO 
C  LBDORI  -  LEFT  BOUNOARY  OENSITY  RATIO  AT  TIME  ZERO 
C  LBOPR  -  LEFT  BOUNOARY  PRESSURE  RATIO 
C  LBOPRI  -  LEFT  BOUNOARY  PRESSURE  RATIO  AT  TIME  ZERO 

C  LBOPRS  -  VALUE  OF  0  DENOTES  CONSTANT  PRESSURE  AT  LEFT  BOUNDARY 

C  ,1  DENOTES  ADJUSTABLE  PRESSURE  AT  THE  LEFT  BOUNOARY 

C  LBOTR  -  LEFT  BOUNOARY  TEMPERATURE  RATIO 
C  LBOTRI  -  LEFT  BOUNOARY  TEMPERATURE  RATIO  AT  TIME  ZERO 
C  LBNDRY  -  DENOTES  LEFT  BOUNOARY  CONDITION,  OPEN  OR  CLOSED 

C  LNODE  -  ARRAY  OF  LEFT  MOST  NOOE  TO  BE  CORRECTED  IN  CORRCT 

C  LHPRES  -  DENOTES  HHICH  SIDE  OF  OIAPHRAGM  HAS  LOH  PRESSURE 

C  N  NUMBER  OF  SPACIAL  NODES  (OOO  NUMBER) 

C  NO  DOUBLE  PRECISION  VALUE  OF  N 

C  NEM**(I)-  STOREO  VALUES  OF  **  FOR  THE  NEXT  TIME  LEVEL 
C  PARRAY  -  ARRAY  OF  PRESSURES  FOR  PLOTTING 

C  PRI  -  INITIAL  PRESSURE  RATIO  ACROSS  THE  SHOCK 

C  PLTCNT  -  COUNTER  FOR  GRAPHICS  ROUTINES 

CO-  ABSOLUTE  FLUID  VELOCITY 
C  QARRAY  -  ARRAY  OF  VELOCITIES  FOR  PLOTTING 
C  QLBD  -  INITIAL  VELOCITY  AT  LEFT  BOUNDARY 
C  QLI  -  INITIAL  VELOCITY  LEFT  OF  THE  DIAPHRAGM 

C  ORBO  -  INITIAL  VELOCITY  AT  RIGHT  BOUNDARY 
C  ORI  -  INITIAL  VELOCITY  RIGHT  OF  THE  DIAPHRAGM 

C  00  Q*A*S  I  EXTENDED  RIEMAFM  VARIABLE ) 

C  RBOOR  -  RIGHT  BOUNDARY  OENSITY  RATIO 

C  RBOORI  -  INITIAL  RIGHT  BOUNOARY  DENSITY  RATIO 

C  RBDPR  -  RIGHT  BOUNDARY  PRESSURE  RATIO 

C  RBOPRI  -  INITIAL  RIGHT  BOUNOARY  PRESSURE  RATIO 

C  RBOPRS  -  VALUE  OF  0  DENOTES  A  CONSTANT  PRESSURE  AT  RIGHT 
C  BOUNOARY,  WHILE  1  IS  FOR  ADJUSTABLE  PRESSURE 

C  RBOTR  -  RIGHT  BOUNDARY  TEMPERATURE  RATIO 
C  RBOTRI  -  INITIAL  RIGHT  BOUNOARY  TEMPERATURE  RATIO 

C  RBNORY  -  DENOTES  RIGHT  BOUNOARY  OPEN  OR  CLOSED 

C  ANODE  -  ARRAY  FOR  RIGHT  MOST  NOOE  TO  BE  CORRECTED  IN  CORRCT 
C  RR  G-A*S  (EXTENDED  RIEMANN  VARIABLE) 

C  S  ENTROPY 

C  SARRAY  -  ARRAY  OF  ENTROPY  FOR  PLOTTING 
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C  SIGMA  -  SPATIAL  LOCATION  OF  DISCONTINUITIES 
C  SIGMA l L ,  J  )  WHERE  L  INCICATES  THE  TYPE  OF 

C  DISCONTINUITY  AND  J  INDICATES  THE 

C  TIME  LEVEL)  1  -  CURRENT  LEVEL 

C  2  -  LEVEL  BEING  CALCULATED 

C  SK  -  INTEGER  THAT  DENOTES  RELATIVE  LOCATION  OF  SHOCK  NEAR 
C  BOUNDARIES 

C  SXIP  -  VARIABLE  WHICH  INDICATES  HOW  MANY  TIME  STEPS  BETWEEN 
C  CALLS  TO  OUTPUT  ROUTINES 

C  T  TIME  SINCE  INITIAL  CONDITIONS 

C  TRI  -  INITIAL  TEMP.  RATIO  ACROSS  THE  SHOCK 

C  VHEAD  -  VELOCITY  OF  THE  HEAD  OF  THE  EXPANSION 

C  WAVE  FOR  THE  EXACT  SOLUTION. 

C  VTAIL  -  VELOCITY  OF  THE  TAIL  OF  THE  EXPANSION 
C  WAVE  FOR  THE  EXACT  SOLUTION. 

C  VCOE  -  VELOCTIY  OF  THE  CONTACT  DISCONTINUITY 
C  FOR  THE  EXACT  SOLUTION. 

C  VS  -  THE  SHOCK  SPEED! POSITIVE  RIGHT,  NEGATIVE  LEFT) 

C  VSE  VELOCITY  OF  THE  SHOCK  FOR  THE  EXACT 

C  SOLUTION. 

C  W  MACH  NO.  RELATIVE  TO  A  STANDING  SHOCK 

C  XARRAY  -  ARRAY  OF  SPATIAL  POSITIONS  FOR  PLOTTING 

C  XEXACT  -  ARRAY  OF  SIX  X  VALUES  FOR  THE  EXACT  SOLUTION. 

C  XINIT  -  INITIAL  POSITION  OF  DISCONTINUITY  FOR 
C  EXACT  SOLUTION  PLOTTING. 

C  X2  LOCATION  OF  NOOE  TO  RIGHT  OF  DISCONT. 

C  ALONG  THE  SPACIAL  AXIS. 

C  Y  <N+l)/2 

C  YEXACT  -  ARRAY  OF  SIX  DENSITY  VALUES  FOR  THE  EXACT  SOLUTION. 

C 

C  ***  OTHER  VARIABLES  ARE  DEFINED  IN  THE 

C  SUBROUTINES  WHERE  THEY  ARE  USED  »*** 

C 

C  ♦♦♦♦♦♦♦♦  ♦♦♦♦♦♦♦♦♦♦ 

c  ♦  ♦ 

C  ♦  PROBLEM  SET-UP  ♦ 


THE  PARTICULAR  PR08LEM  FOR  THIS  VERSION  IS: 

SHOCK  TU8E ,  SINGLE  CENTERED  OIAPHRAGM  WITH 
HIGH  PRESSURE  SIDE  TO  THE  RIGHT. 


BOUNOARY  CONDITIONS  -  LEFT  END  CLOSED, RIGHT  END  OPEN 


VARIABLE  DECLARATIONS 


DIMENSION  1 2 ( 4 )  ,X2( 4  )  .XEXACT ( 6  )  , YEXACT ( 6  ) , LNODE ( 4  )  ,RNOOE ( 4 1 , 
C  SIGMA! 4,2  I 

USER  INPUT  REQUIRED  HERE  ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 


SET  THE  DIMENSIONS  EQUAL  TO  N 


DIMENSION  A( 101 ) ,Q< 101 ) ,QQ( 101 )  ,RR( 101 )  ,SI 101 ) , 
C  NEWRR! 101 ), NEWS! 101), NEWQQ! 101 ), 

C  PARRAYI 101 J ,D ARRAY! 101 )  ,SARRAY( 101 ) , 

C  QARRAY! 101), XARRAY!  101) 


INTEGER  I , J ,N , JSTOP ,Y .GRAPHS .COUNT .PLTCNT ,B0RY ,SK .RBNDRY , 

C  SKIP  ,12  ,  LWPRES .HALT , LNODE .RNODE , LBNORY , LBOPRS .RBOPRS 

DOUBLE  PRECISION  TRI ,PRI ,QLI ,QRI ,0RI ,G,G1 ,G2 , SIGMA , EE , NEWQQ, 

C  0ELT,H,N0,X2,AR,W,0Q,VS,T,A,Q,QQ,RR,S,  NEWRR,  NEWS, 

C  LBOPRI , LBOTRI , LBDORI , LBDPR . LBOOR . LBOTR ,QLBO , 

C  RRA ,QQA  > AA ,5A ,QA .RRB ,QQB . AS ,SB ,QB , 

C  RBOPRI  ,RB0TRI  .RSODRI  .RBOPR , RBOOR .RBOTR ,QRB0 , 

C  QCS.VTEW.VHEW 

REAL  VTAIL , VCOE ,VSE ,DLSH ,DLC0 , XINIT , VHEAD , 
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C  XEXACT , YEXACT *P ARRAY ,0 ARRAY .SARRAY ,9 ARRAY .XARRAY 

comoN  AR.DG.VS.W 
C 

C  -  ENTER  THE  APPROPRIATE  VALUES  BE LON  - 

C 

N»  101 
GRAPHS*! 


C 

C  -  FOR  GRAPH  *  1  OR  2  MUST  ENTER  CHANGES  IN  SUBROUTINE 

C  -  BORDER  AND  SUBROUTINE  EXACT  - 

C  -  LINES  "FIRST  ORDER  N  «  ???" 

C  -  "DENSITY  RATIO  *  ???  TEMP  RATIO  *  ???" 

C  -  "PRESSURE  RATIO  *  ???" 

C 


SKIP* 18 
JSTOP*101 
TRI*1.0000 
PRI*5. 0000 
ORI*S. 0000 
QLI*0 . 0000 
QRI*0. 0000 
LBOPRI  *  1.000 
LB0TR1  *  1.000 
LBOORI  *  1.000 
RBOORI  *  1.000 
RBOPRI  «  4.000/5.000 
RBOTRI  a  1.000 
Gsl.4000 
EE-0.1D-8 
C 


C  - -  DENOTE  LON  PRESSURE  SIDE  BY  SETTING  - 

c  -  LWPRES  a  2  IF  LON  PRESSURE  ON  RIGHT 

C  -  LWPRES  a  3  if  LON  PRESSURE  ON  LEFT 

C 

LWPRES  a  3 
C 

c  -  SET  BOUNDARY  CONDITIONS  BY  SPECIFING  OPEN  OR  CLOSED  - 

C  -  LBNDRY :  OPEN  *  0,  CLOSED  a  1  (FOR  LEFT  BOUNDARY)  - 

C  -  IF  OPEN  SPECIFY  IF  PRESSURE  IS  TO  BE  MAINTAINED  AT  LBOPRI  - 

C  -  OR  IF  IT  CAN  AO JUST  TO  PREVENT  ANY  WAVES  FORMING  AT  THE  BOUNDARY  - 

c  -  LBOPRS:  CONSTANT  *  0,  ADJUSTABLE  *  1  - 

C  -  DO  THE  SAME  FOR  RIGHT  BOUNOARY \  RBNORY , RBOPRS  - 


C 

LBNDRY  a  1 
LBOPRS  a  1 
RBNORY  a  0 
RBOPRS  *  0 

C  -  EXACT  SOLUTION  VALUES  - - 

XINITa0.50 
VHEAD*  1.0 
VTAIL*  0.510557 
VCDE»-. 574487 
VSEa-1. 402546 
0LC0*2. 715115 
DLSH*1. 64544 
SIGMA ( 1 ,2  )*o . 50000001000 
SIGMA! 2 , 2  )*o . 50000001D00 
SIGMA! 5,2 1=0.50000001000 
SIGMA! 4 ,2  >=0 . 50000001000 
C 

C  END  OF  USER  INPUT  AREA  ♦♦♦♦♦♦♦♦ 

C 

SK  ■  0 
T=0 . 0000 
00  10  1=1,4 
12!  I  >=0 
LNOOEII)  *  0 
RHODE  ( I )  *  0 
X2(  I)=0.0000 


uuuuuu  ouu  o  u  ooo 


HIUIUIUIUIIWIHIUIIlP'JilUMUi 


SIGMAd.l  1*0.0000 

10  CONTINUE 

N0*0BLE(N1 

H*1 . D00/1 N0-1.000 ) 

00  11  1*1 »N 

A(  I  1*0.0000 
Q( 1 1*0.0000 
NEHQOI I >*0.0000 
NEHRRI I ) *0.0000 
NEHS( 1 1*0.0000 
QARRAYI 1 1*0.0 
P ARRAY! 1 1*0.0 
D ARRAY (1 1=0.0 
SARRAYII  1=0.0 

XARRAY <  1 1=FL0ATI  1-1 1»SNGL<  H 1 

11  CONTINUE 

0ELT*2.0000 

- LOAD  INITIAL  REIMAN  VALUES  INTO  NODE  LOCATIONS .FIRST  FROM  NODE— - 

- 1  TO  MIOPOINT  (Yl,  ANO  THEN  FROM  Y  TO  N.  NOTE  IF  SHOCK  DOES  NOT— 

- START  AT  MIDPOINT  THEN  Y  SHOULD  BE  SET  TO  NODE  WHERE  SHOCK  - 

- INITIALLY  IS - 

AR* 1.0000 

DQ*0.0000 

H*1.0000 

VS*0.0000 

QCS*0.000 

VHEM*0.000 

VTEW*0 . 000 

01*1. 000/1  G-l. 000) 

62*2. 000/1 G-1.D001 

-  FLOW  RIGHT  - 

IF<  LWPRES.EQ. 2  1  THEN 
LBOOR  *  DRI  »  LBDDRI 
LBOPR  *  PRI  *  LBOPRI 
LBOTR  *  TRI  *  LBOTRI 
QLBD  *  QLI 
RBOOR  *  RBOORI 
RBOPR  *  RBDPRI 
RBOTR  *  RBOTRI 
QRBO  *  OR I 

Y  *  (  N*l  1/2 
Y  *  38 

DO  12  1*1. Y 

SII  »=G2-t  Gl/G  1*OLOG<  PRI/I I  DRI  1**G  1 1 
QQ( I  )=QH*0SQRT( TRI  )*Sl 1 1 
RR(  I  l=QLI-OSQRT(  TRI  )»S(  I ) 

12  CONTINUE 

Y«tN*3  1/2 

Y*39 

00  II  I*Y,N 
SI  I  l*GZ 
6QII  )*QRI«Sd  1 
RRI 1 1*QRI-SI  I  1 

13  CONTINUE 
ELSE 

-  FLOW  LEFT  - 

L800R  *  LBDDRI 
LBOPR  *  LBOPRI 
LBOTR  *  LBOTRI 
QLBD  *  QLI 

RBOOR  *  RBOORI  *  ORI 
RBOPR  *  RBOPRI  »  PRI 
RBOTR  *  RBOTRI  «  TRI 
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c 


QRBO  •  QRI 

Y  «  IN-1J/2 

Y  ■  13 

00  14  1*1 ,Y 
SI  I  )*G2 

QQI I  )=QLI*S(  I  ) 

RRI I  )*QLI-S( I ) 

14  CONTINUE 

Y  «  INUJ/2 

Y  *  14 

DO  IS  I*Y,N 

SI  I  )=G2-I  Gl/G  )»0L0G1  PRI/I I  DRI  )**G ) ) 

QQI I  )=QRI+OSQRT( TRI  l*S! I  ) 

RRI I 1=QRI-0SQRT l  TRI )»SI  I ) 

15  CONTINUE 
END  IE 

-  SET  UP  GRAPHICS  PLOTS  OF  VARIABLES  - 

IF  I  GRAPHS. GT.O)  THEN 
CALL  COMPRS 
CALL  TEK618 

CALL  HWROTI 'AUTO' ) 

CALL  HWSCALI 'SCREEN* ) 

IF  I  GRAPHS. EQ. 2  )  THEN 
CALL  PAGE! 11.0,8.5) 

ELSE 

CALL  PAGE 1 8. 5 >11.0  ) 

END  IF 

IF  I  GRAPHS. EQ.l)  THEN 
CALL  BORDER! JSTOP ) 

END  IF 
ENO  IF 
HALT  *  0 
J=1 

C0UNT*1 

IF  I  GRAPHS. EQ.l)  THEN 

CALL  P LOT I J , JSTOP , N , QQ , RR , S , H ,XARRAY , PARRA Y , 

RDARRAY .QARRAY >SARRAY ,G , G1 ,G2 ) 

END  IF 

-  BURST  DIAPHRAGM  - 

CALL  TIME(N,Qq,RR,S,DELT,Hl 

CALL  OBURSTI N ,H ,QQ ,RR ,S,G,G1 ,G2  ,DELT  tI2  >X2 ,H, AR ,DQ,VS,LNPRES  , 

C  SIGMA, A, Q I 

IF  I  GRAPHS. EQ.O)  THEN 

CALL  LISTIN, SIGMA, QQ,RR,S,G,G1,G2,J,T,DELT, VS, QCS.VTEW.VHEH, 
C  XINIT ,VHEAD ,VTAIL ,VCDE ,VSE ) 

END  IF 

-  BEGIN  CALCULATION  FOR  JUMP  TO  NEXT  TIME  AND  CONTINUE  - 


UNTIL  EITHER  JSTOP  REACHED  OR  SHOCK  MEETS  CONTACT  SURFACE 


16  IF  I J.EQ. JSTOP )  GOTO  18 
PLTCNTsJ/SKIP 


CALL  TIME! N,QQ,RR,S, CELT ,H I 


i 


CALL  TRAKI N, SIGMA ,H,QQ,RR ,S,G»G1 ,G2 ,0ELT,I2,X2,N,AR,DQ,VS,J, 
C  LHPRES ,QCS,VHEH ,VTEW I 


CALL  SWEEP! N, H, SIGMA, Qq,RR,S,DELT*EE,Q, A, NEWQQ.NEWRR, NEWS, 12, G2, 
C  J .LNODE ,RNODE .HALT .LBNORY .LBDPRS . LBOPR .LBDTR , LBOOR . 

C  QLBO.G.Gl.RBNORY .RBOPRS .RBDPR .SBOTR ,SBDOR .QRBO ,3DRY , 

C  SK  I 


IFI LNODE 1 4 ) . LT . J )  THEN 


i 

I 


1 


»*■  Y.v  JvV. V-  XL  i£i. V 


LNOOEUI  ■  0 
RNOOEIl)  •  0 
GO  TO  17 

END  IP 

XFl RNOOEt  4 ). LT . LNOOEI 4) )  THEN 
RNOOEIl)  ■  0 

ENO  IF 
C 

CALL  CORRCTI  LNODE  .RNODE  >N,SIGMA,H ,Q9,RR,S ,G ,G1 ,G2 ,12 >X2 ,H,AR ,DQ, 

C  VS, A, 9) 

C 

17  IF  ( SIGMA ( 1 ,2 ) .GE .1.000 )  THEN 
IF( RBNORY . EG. 0 )  THEN 

IFISIGMAll,2).NE.3.D00)  THEN 
SK  »  2 

CALL  B0NDRYI9IN),9IN'1),QRBD,AIN>,AIN-1),QQIN),99IN-1), 
C  RR(N),RR(N-1), SIN), SIN-1  ),H,EE,DELT, 

C  RBNORY iRBDPRS ,RBOPR ,RBOOR ,RBOTR ,J ,NEHQQ( N ) , 

C  NEHRRI  NI,NEHS(N),G,G1>G2  >HALT  ,  BDRY  ,SK  ) 

ENO  IF 

ELSE 

CALL  EXTRAPI RR( N-l ) ,RRI N-2  )  ,GGI N-l ) ,QGI N-2 ) ,SI N-l ) , 

C  SIN-2  ),H,H,RRA,QQA,SA,AA,QA  ) 

CALL  SRFLCTI QQA.RRA ,SA , SIGMA , VS , CELT , LHPRES, 

C  RR(N),QG(N),S(N),QIN),AIN),G,G1,G2) 

ENO  IF 

ENO  IF 
C 

IF  ( SIGMA I  1 , 2 ) . LE . 0 . 000 )  THEN 
IF! LBNORY.EQ.O I  THEN 

IF( SIGMA! 1,2  I. NE. -2.000  )  THEN 

BORY  *  3 

SK*2 

CALL  BONORYI 91 1 1  ,Qt  2 1  ,QLBD  ,AI  1 ) , Al  2  )  ,091  1 )  ,091  2  ) , 

C  RRI 1 )  ,RRI  2  )  ,S<  1 )  ,S(  2  )  ,H  ,EE  ,OELT , 

C  LBNORY , LBOPRS , LBOPR , LBOOR , LBOTR , J ,NEH99< II, 

C  NEHRRI 1 ) ,NEWS( 1 ) ,G ,G1 ,G2 .HALT ,BDRY ,SK  ) 

ENO  IF 

ELSE 

CALL  EXTRAPI  RRI  2  )  ,RR{  3  )  ,091  2  )  ,091  3  )  ,S(  2  ) , 

C  SI  3  )  ,H  ,H  ,RRB  ,008  ,S8  ,AB  ,QB  I 

CALL  SRFLCTI QGB ,RR8 ,SB .SIGMA , VS ,DELT .LHPRES , 

C  RRI  ll.QQI  1  ),SI  11,01  1  >,A(  1  I,G,G1,G2  ) 

ENO  IF 

END  IF 

T=T+DELT 

C 

C - OUTPUT  DATA - 

C 

IF  I  I  COUNT . EG . PLTCNT»SKIP  ) . AND . I  GRAPHS . EQ . 1  )  ) 

C  THEN 

C  IFI I J.GT.55 ) )  THEN 

CALL  PLOT! J ,JSTOP ,N,Q9»RR »S»H,XARRAY .PARRAY , 

CDARRAY ,QARRAY .SARRAY ,G ,G1 ,G2  I 
C  ENO  IF 

END  IF 

IF  1 1  COUNT . EO . PLTCNT«SKIP J . AND . I  GRAPHS . EG . 0  I) 

C  THEN 

CALL  LIST) N, SIGMA, 09, RR,S,G,G1,G2,J,T,DELT, VS, 9CS,VTEH,VHEH, 

C  XINIT ,VHEAD .VTAIL ,VCDE ,VSE  ) 

ENO  IF 

IF  I  I  COUNT  . EO. PLTCNT*SKIP ) .AND. I  GRAPHS. EQ. 2 ) ) 

C  THEN 

IF  IJ.GT.50)  THEN 

CALL  EXACTI N .XINIT ,T .VHEAD , VTAIL ,VCDE ,VSE ,DLCD,OLSH ,99 ,RR ,S ,H , 
CXARRAY  ,0 ARRAY  ,G,G1,G2,0RI  I 

ENO  IF 

ENO  IF 
C 
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onoooonoooooo 


inHALT.BQ.il  GO  TO  18 
J*J*1 

COUNT "COUNT +1 
GO  TO  16 


CALL  OONEPL 
END 


SUBROUTINE  LISTI N, SIGMA ,QG,RR ,S,G,G1 ,G2 , J,T ,DELT ,VS .GCS.VTEH, 
C  VHEH.XINIT ,VHE AD ,VTAIL ,VCDE ,VSE  ) 


♦  TABULAR  RESULTS  SUBROUTINE  ♦ 


♦  * 


VARIABLE  DEFINITIONS 


DENS  -  DENSITY 
PRESS  -  PRESSURE 
TEMP  -  TEMPERATURE 

INTEGER  I,J,N,L 

DIMENSION  SIGMA(4,2I,QQINI,RRIN1,SIN1 

DOUBLE  PRECISION  SIGMA  ,QQ,RR,S, PRESS,  VTEH.VHEN.QCS. TE»«E , 

C  TEMP, DENS, G.G1.G2.Q.T.DELT, VS, HEHXE, 

C  XINIT ,VHEAD .VTAIL ,VC0E , VSE ,SKXE ,CSXE 

HRITEI  9,*  I  'TIME  LEVEL1, „V  ELAPSED  TIME  IS'.T 

HRITEI  9>»  I  'TIME  STEP  IS', CELT,'  SHOCK  VELOCITY  IS’, VS 

NRITEI 9,* )  'CONTACT  SURFACE  VELOCITY  IS'.QCS 

HRITEI 9, »  I  'HEAD  EXPANSION  HAVE  VELOCITY  IS’ ,VHEH 

HRITEI 9, *  I  'TAIL  EXPANSION  HAVE  VELOCITY  IS'.VTEH 

HRITEI  9,*;  '  ’ 

HRITEI 9, • I  •  NOOE  VELOCITY  DENSITY 

CPRESSURE ' 

HRITEI 9.*  I  *  ’ 

00  A1  1*1, N 

TEMP*!  OQI I  l-RRII  I  1*1  QQI I  l-RRI  1 1  l/<  4.D00»SII  )*SI  1 1 1 
DENS*l  11.000/TEMP  )»OEXP( G*(  1 . 000-G  1*1  St  I  )-G2  I I  1**1  -G1 1 
PRESS*TEMP»OENS 
Q*IQQI  I  )*RR(  I  )  1/2. 0000 
HRITE  I  9,65  I  I ,Q,OENS, PRESS 
AS  FORMAT  I 4X.I2 ,7X,F12 . A ,7X,F12 .A ,7X,F12.A I 
A1  CONTINUE 

HRITEI*, »>  '  * 

HRITEI 9,*)  '  ' 

HRITEI 9,*  I  '  NOOE  GO  RR 

CMOOIFIED  S' 

HRITEI 9,* I  ' 

DO  A2  1*1. N 

HRITE  I 9,AA  I  I ,QQ< 1 1 ,RRI 1 1  ,S<  1 1 
AA  FORMAT  I 4X.I2 ,7X ,F12 . A .7X.F12 . A ,7X ,P12 . A  I 
A2  CONTINUE 

HRITEI 9,*  I  ' 

HRITEI 9,*  I  '  • 

HRITEI 9, » I  '  DISCONTINUITY  LOCATIONS  AT  TIME  LEVEL', J 

HRITEI  9.* I  •  • 

HRITEI 9, »  I  '  TYPE  LOCATION* 

DO  A3  1*1,4 

HRITEI 9,*  I  L , '  1 .SIGMA I L ,2  I 

43  CONTINUE 

HRITE'*.*'  '  ’ 

HRITEI  9, »  I  ’ - - - - - 

C----- . — . . . . 

miTE<9,»l  '  ' 
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OOUUUO<JUU<J<J  o  uouuuuuuouuuuuuu 


IFIJ.E4.1)  THEN 

MRXTEI 10,* )  *  TIME  SHOCK  CONTACT  HE AO 

C  TAIL* 

MRXTEI 10,*)  '  * 

I  NO  If 

MRXTEI 10,471  T , SIGMA! 1,1)  , SIGMA! 2,1) ,SXGMA( 5 ,1 ) .SICMAI  4,1) 

47  FORMAT  I1X,F12.4>1X»F12.4,1X,F12.4,1X,F12.4,1X»F12.4) 
XFIJ.E4.il  THEN 
NRITEIO,*)  '  EXACT  VALUES' 

MRXTEI*.* )  *  TIME  SHOCK  CONTACT  HE AO 

C  TAIL' 

MRITEl  *.» I  '  ' 

END  IF 

SKXE«XINIT«VS£*T 
CSXE*XXNIT ♦VCOE*T 
HEHXE  «XINI  T  ♦' VHE  AO*T 
TEMXI *XINIT*VTAIL»T 
MRXTEI  0,4*  I  T  ,SKXE  >CSXE  .HEMXE  »TEW<E 
40  FORMAT  1 1X,F12.4,1X,F12.4,1X,F12.4,1X,F12.4,1X,F12.4  ) 

RETURN 

ENO 

C 

SUBROUTINE  TIME! N.44 ,RR,S, CELT ,H ) 


a  CALCULATE  TIME  STEP  SUBROUTINE  a 


- NEM  VARIABLE  DEFINITIONS  - 

THIN  -  RLM4ING  VALUE  OF  THE  MINIMUM  TIME  STEP 
INTEGER  N.I 

DIMENSION  QOIN).RR(N).SIN) 

OOUBLE  PRECISION  H,A,4Q,RR,S,0ELT,TMIN,4 

TMIN*2.0000 

DO  21  1*1 ,N 

A*IOQI  I  l-RRI  I ) )/(  Z.DOOaSI  I ) ) 

QXQQII  )*RRI  I  I  1/2  .  ODOO 
OELT*H/t DABSI OABSI Q  l*A  I  I 
IF  IOELT.LT. TMIN I  THEN 
TMIN*OELT 
EM)  IF 
21  CONTINUE 

0ELT*0 . TTOOO»TMIN 

RETURN 

ENO 

SUBROUTINE  TRAK IN, SIGMA ,H ,QQ,RR ,S ,G ,G1 ,G2 ,DELT ,12 ,X2 ,M,AR ,04, VS 
OJ,LMPRES,QCS,VHEM,VTEM ) 


a 


a 


a  DISCONTINUITY  TRACKING  SUBROUTINE  • 

a  a 


VARIABLE  DEFINITIONS 


CSOXR 

CSRFtl 

OR 

0REM4 

PR 


■  CONTACT  SURFACE  DIRECTION,  2  TO  THE  RIGHT,!  TO  THE  LEFT 
'  RXEMANN  VARIABLE  CHANGE  ACROSS  A  CONTACT  SURFACE 

THE  RATIO  OF  THE  DENSITY  ACROSS  A 
SMOCK,  3/AI  RIGHT  I  ,B/AI  LEFT  I 

■  DUMMY  VARIABLE 

THE  RATIO  OF  THE  PRESSURE  ACROSS  A 
SHOCK,  A/SI  RIGHT  )  ,0/AI LEFT  ) 


nooo  non  noon  ono  ononoooooooo 


HREIMN  -  THE  MEASURED  JUMP  IN  09  ACROSS  THE  SHOCK, 

FROM  A  TO  B. 

EREIttt  -  THE  JUMP  IN  qQ  ACROSS  THE  SHOCK  CALCULATED  ANALYTICALLY 
AS  A  FUNCTION  OF  H. 

SAP  -  ENTROPY  TO  THE  LEFT  OF  THE  SHOCK  FOR  FLOWS  RIGHT 

OR  ENTROPY  TO  THE  RIGHT  OF  THE  C.S.  FOR  FLOMS  LEFT 
S8P  -  ENTROPY  TO  THE  RIGHT  OF  THE  C.S.  FOR  FLOWS  RIGHT 

OR  ENTROPY  TO  THE  LEFT  OF  THE  SHOCK  FOR  FLOWS  LEFT 
SHKDIR  -  SHOCK  DIRECTION  OF  TRAVEL, 3  TO  THE  LEFT, 2  TO  THE  RIGHT 
TS  -  TIME  FOR  SHOCK  TO  TRAVEL  ONE  INTERVAL 
X  -  DISTANCE  FROM  LEFT  BOUNDARY  TO  NODE 

INTEGER  N,I»Y>I2»L,J , SHKDIR »CSOIR >LWPRES 
DIMENSION  SIGMA! 4, 2 1  ,X2!  4  )  ,RR(  N  )  .SQ!  N  )  ,SI  N  )  ,12!  4  ) 

DOUBLE  PRECISION  SIGMA ,X2 ,X,H ,AB ,SA , SB ,AA ,QA ,QB,QQA , GOB ,RRA ,RRB, 
C  RR,QQ,S,TS,M,DQ,AR,PR,G,G1,G2,VS,0ELT,CSRMN, 

C  G,MREIMN .OREMN .EREIMN ,WW ,SAI ,SA2 ,SAP ,SBP 

C  AH,QH,VHEH,VTEH,TIME.QCS 

♦♦444444  LOCATING  THE  UPSTREAM  NOOE  ♦♦♦♦♦♦♦♦♦♦♦ 

DO  10  L«l,4 

SIGMA! L,1  )*SIGMA( L ,2  ) 

Y«0 

XsO.OOO 

I«1 

11  IF  !  . NOT. (Y. EG. OH  GOTO  10 

IF  I  SIGMA! L,1 1 . LT .X )  THEN 
X2ILI«X 
I2(L)»I 
Y=1 

END  IF 
X*X+H 
I«I*1 
GOTO  11 
10  CONTINUE 

-  IF  SHOCK  HAS  LEFT  AN  OPEN  BOUNDARY  OUT  OF  THE  TUBE  THEN  SET  - 

-  SHOCK  TO  NEUTRAL  - 

IF  112(1). GT.N)  THEN 

SIGMA! 1,1)  3  z.ooo 
SIGMA! 1,2)  a  3.000 
W  s  1.000 
VS  =  1.D00 

OQ  3  0.000 

AR  «  0.000 
PR  a  1.000 
DR  3  1.000 
GO  TO  15Q 

ELSE  IFII2ID.LT. 2)  THEN 
SIGMA! 1,1)  a-1.000 
SIGMA! 1 ,2  )  s-2.000 

H  a  1.000 

VS  a-l.DOO 
OQ  a  0.000 
AR  a  0.000 
PR  »  1.000 
OR  ■  1.000 
GO  TO  ISO 

END  IF 

4444444444444  DETERMINING  SHOCK  SPEED  44444444 

IF!  ( J.EQ.  1 1.OR.  ( 12!  D.E0.2  I. OR.  ( 12!  1  I.EQ.N)  )  THEN 

- AT  TIME  ZERO  OR  BOUNOARYS  OETERMINE  CORRECT  SHOCK  DIRECTION - 

- SHKDIR  a  J  is  A  SHOCK  HEADED  LEFT,  AND  SHKDIR  ■  2  IS  SHOCK - 

- TRAVELING  RIGHT - 
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conn  o  noon  noon  ooonn 


C 


IF(LNPRES.EQ.S)  THEN 
SHKDIR  ■  3 
IF  ( J.EQ.l)  THEN 

X2( X  )  ■  X2(l>  -  H 
12(1)  ■  12111  -  1 
X2(  2 )  •  X2(  2 1  -  H 
12(21  ■  12(21  -  1 

END  IF 
60  TO  20 

ELSE 

SHKDIR  «  2 

END  IF 
SO  TO  20 

—IF  SHOCK  ANO  CONTACT  SURFACE  ARE  NOT  WITHIN  THE  SAME - 

— INTERVAL  THEN  NO  CORRECTIONS  ARE  NEEDED  IN  CALCULATING - 

— THE  REIMAN  VARIABLE  JUMP  ACROSS  THE  SHOCK - 

ELSE  IF  (12(11. NE. 12(21)  THEN 

IF(  ( SHKDIR.  EQ.  31.  AND.  (SIGMA!  1 ,1 1.  EQ.  I  X2(  1  )-H  111  THEN 
X2(  1 1  ■  X2(  1 1  -  H 
12(11  *  12(11  -  1 

END  IF 
GO  TO  20 

— IF  SHOCK  ANO  CONTACT  SURFACE  ARE  WITHIN  THE  SAME  INTERVAL - 

— THEN  CORRECTIONS  ARE  REQUIRED  TO  DETERMINE  SHOCK  STRENGTH - 

ELSE  IF( SHKDIR. EQ. 31  THEN 

— SHOCK  LOCATION  RELATIVE  TO  THE  CONTACT  SURFACE  FOR  A  SHOCK - 

- HEADED  TO  THE  LEFT  DETERMINES  THE  CORRECT  VALUES  FOR  H  AND  VS— 

IFf  SIGMA! 1 > 1 1 .GT. SIGMA!  2,11)  THEN 
GO  TO  111 

ELSE  IF (SIGMA 1 1,1 1. EG. 1X2 ( 1 )-H 1 1  THEN 
X2( 1 1  «  X2(ll  -  H 
12(11  «  12(11  -  l 
X2(  2 1  ■  X2(21  -  H 
12(2)  *  12(2)  -  1 

END  IF 

IS  RRA=RR(  12(1 1-1) 

RRB=RR(  12(1)1 
QQA=qq<  12(1  i-i) 
qqs=qq< 12(11) 

SA=S<  12(1 1-1 1 

QA  *f QQA+RRA  1/2 . DOO 
QB  =  <  QQB+RRB  1/2.000 
AA  s<  OQA-RRA  )/( 2 . DOO*SA  1 
Oq  »<QB-QAI/AA 

W  *  OSQRTI ( DQ»»2 )*( 0.36D00  1+1 . DOO 1  -  (DQ*0.6000) 

Oq  *<  -1.000  l»OQ 

GO  TO  110 

- FOR  SHOCK  HEADED  RIGHT  THE  SAME  PROCEDURE  FOR  CORRECTIONS  ARE - 

- FOLLOWEO - 

ELSE 

ELSE 
17 
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IFISIGMAI 1,1). LT.SIGMAI 2,1)1  THEN 
GO  TO  111 

RRA*RR(I2<1)-1) 

RRB=RR( 12(1)1 
QQA*QQ(  12(1 1-1) 

QQB=QQ(  I2(  1 1 ) 

SB-31 12(1 1 1 
qA  *(QQA*RRA  1/2. DOO 
QB  ■(qqe»RRB)/2.D00 


n  n  r>  r>  nnnnn  o  n  o  o  o  n 


AB  »( QQB-RRB )/( Z . D00*SB ) 
tXJ  a(QA-QB)/AB 

N  a  DSQRTI ( 0Q**2 )*( 0.36000 1+1 .DOO )  ♦  (DQ«0.6000> 
60  TO  110 

END  if 

- KITH  NO  SHOCK/CONTACT  SURFACE  INTERACTION  THE  JUMP  IN  REIMAN - 

- VARIABLES  ARE  DETERMINED  WITHOUT  INTERPOLATION  OVER  THE  -  - 

- INTERVAL.  MREIMN  IS  THE  MEASURED  JUMP.  EREIMN  IS  THE  ANALYTICAL-- 

- VALUE - 


20  RRAaRRl  121 11-1) 

RRBsRRI I2( 1 1 ) 

QQAsQQI  121 1  >-l) 

QQBsOOC I2( 1 ) ) 

SA»SII2il)-l) 

SB»S(I2(in 

21  AB>( QqB-RRB )/( 2.D00*SB ) 

AA*<  OQA-RRA  )/(  2 . 000*SA  I 
IF( SHKDIR . ED. 3 )  THEN 

MREIMN  >  ( RRB-RRA  )/AA 
DREMN  >  DABSI  MREIFN) 

ELSE 

MREIMN  >  (QQA-QQBI/AB 
DREMN  «  MREIMN 

END  IF 

- ITERATE  FOR  PROPER  VALUE  OF  N  USING  THE  QUADRATIC  FIT  OF  THE- 

- REIMAN  VARIABLE  CHANGE  WITH  W  CURVE.  NOTE  LEFT  MOVING  SHOCKS 

- ARE  USED  IN  THESE  EQUATIONS  SINCE  RRB-RRA/ AAa-( QQA-QQB/AB  ) - 

100  HH»  ( 3. 0396408D01-H0REMN+2. 7574D00  )/0. 28633 7D00  ) ) 
H*5.513294000-0SQRT( Hh 1 
0Q«2 . 000*1 HMW-1 . DOO  l/( H*(  G+l . DOO  )  ) 

ARaOSQRTI  2 .  D00*(  G-l .  000  l*(  1 .  DOO+I  I  G-l . DOO  )*M*W/2. 000  I  )* 

C  I G*G2*H*W-1 . 000  I )/( I  G+l. 000  t*H ) 

PR*( 2.000*G/( G+l . DOO ) )»H*H-( ( G-l . DO0 1/1 G+l . DOO  ) ) 

DR*(  <  G-l .  000  l*H*H+2 .  DOO  I/I  ( G+l .  000  l*H*W  I 
EREIMN=OQ+( AR-l.DOO  )*G2-( AR*G1/G  )*OLOGI PR*( DR**G )  I 
IF  I DABSI EREIMN-OABSI MREIMN 1 ) . LT .0. ID-5 )  GO  TO  110 
DREMN  a  (OABSl MREIMN  I  -  EREIMN)  *  DREMN 
GOTO  100 

- SHOCK  VELOCITY  OEPENDS  ON  DIRECTION  SHOCK  IS  TRAVELING  - 

- LEFT  IS  <  0,  ANO  RIGHT  IS  >  0  - 


110  IF  IJ.EQ.l)  THEN 

TIME  «  0.000 

SA1  a  ( Gl/G  )*0L0G! ( 2 . DOO*G*<  W**2 1-G+l .000 )/( G+l . 000 ) I 
SA2  a  G1*0L0G<  (  (G-l. 000  )*(W**2  )+2  I/I  (G+l. 000  )*(W**2  )  I  I 
IF! SHKOIR. EQ. 2  I  THEN 

SAP  a  SB  -  SA1  -  SAX 
SBP  a  SAP 

ELSE 

SBP  «  SA  -  SA1  -  SAX 
SAP  a  SBP 

END  IF 

103  IF! SHKOIR. EQ. 2  I  THEN 

CSRFN  «l  I DEXPI I SBP-SA  )/G2  1 1*1  SA  l-l  SBP )  >»AR 

ELSE 

CSRMN  ■(  I  DEXPI  I SAP-SB  I/G2  )  )*(  SB  )-( SAP  )  )«AR 

ENO  IF 

IF!  SHKOIR  .EQ.  3  )  THEN 

MREIMN  »(l  RRB-RRA  l/AAI+CSRMN 
DREMN  a  DABSI  MREIMN I 

ELSE 

MREIMN  »l  (QQA-QQB  l/AB  l-CSRff) 

DREW  •  MREIMN 
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W*(3.  0396408001-1 1  DREW *2 . 7574000  )/0 . 286337000 ) 1 
H»S . 51S294000-0SQRT I HH ) 

DQ«2 . 000*1 H*H-1 . 000 )/( H*1  6*1 . 000 ) 1 

AR«DSQRT( 2 . 000*1 6-1 . 000  1*1 1 . 000*1 1 6-1 . 000  )*W*H/2 . 000 ) >* 
C  ( G*62*N*H-1 . 000  )  1/1 16*1 . 000  )*W) 

PR*( 2 . 000*6/1 6*1 . 000 )  )*H*H-( 1 6-1 . 000 1/1 6*1 . D00 ) ) 

DR«( 1 6-1 . 000  )*H*H*2 . 000  »/l  1 6*  1 . 000  >*H*H 1 

EREIW*0Q*( AR-1 .000 )*G2-( AR*G1/G )*OLOG( PR*( 0R»*6 ))  . 

IF  IOABSlEREIMN-OA8S(MREIWn.LT. 0.10-5)  GO  TO  102 
DREW  *  1  DABS!  MREIW )  -  ERE1W)  *  OREW 
GOTO  101 

102  SA1  *  < Gl/G )*DL0GI I 2.000*G*I H**2  1-6*1.000 )/( 6*1.000 1 1 

SA2  >  G1*OLOG(  1(6-1.000  )*(W**2  )*2  )/!  (6*1.000  1*1  M**2  1) ) 
IF( SHKOIR.EQ. 2 )  THEN 

SAP  ■  SB  -  SA1  -  SA2 

ELSE 

S8P  ■  SA  -  SA1  -  SA2 

END  IF 

IF  (OABSISAP-SBPl.LT. 0.10-5)  GO  TO  105 
IF( SHKOIR. EQ. 2 )  THEN 

S8P  «  1 SAP-S8P  )  *  S8P 

ELSE 

SAP  ■  ( S8P-SAP  )  *  SAP 

END  IF 
GO  TO  103 

END  IF 

105  IFI SHKOIR.EQ. 3)  THEN 

OQ  *  (-1.0001*06 

VS  «  ( I RRA*QQA  1*0 . 5000 1  -  1 H»AA 1 

ELSE 

VS*( QQ8*RRB  1*0 .5000 *M*AB 

ENO  IF 

TS*H/DABS( VS 1 
Ill  IF  (TS.LT.OELTI  THEN 
0ELT*0. 99D00*TS 

ENO  IF 

SIGMA! 1,21  *VS*OE LT+SIGMAI 1,11 
C 

C  ♦♦♦♦♦♦OETERMINE  CONTACT  SURFACE  SPEEO********* 

C 

IF(J.EQ.l)  THEN 

CSOIR  *  SHKOIR 

ENO  IF 
C 

c  -  CONTACT  SURFACE  TRAVELING  RIGHT  - 

C 

150  IFI CSOIR. EQ. 2)  THEN 
C 

c  - CONTACT  SURFACE  MOVING  RIGHT,  CHECK  FOR  SHOCK  IN  INTERVAL - 

C  - AND  CALCULATE  SPEEO  OF  CONTACT  SURFACE  AS  APPROPRIATE . 

C 


IFIJ.EQ.il  THEN 

QB  *( QQB*RRB  1/2.000 
QA  *  OQ*AB  *  QB 

CALL  SKJUMPI AB,QB,S8,AR,0Q,VS,G,61 ,H,AA ,QA,SA 1 
Q  ■  QA 
GO  TO  200 


ENO  IF 

IFI  X2(  2  1  .EQ.X2I 1 1 1  THEN 

IFI SIGMA! 2 , 1  I . LE .SIGMA! 1,1)1  THEN 

Q  ■  I QQI 12(21-1)  *  RR( 12121-11)  /  2.000 

ELSE 

Q  *( QQI  1212)1  *  RR(  1212)1)  /  2.000 

ENO  IF 

ELSE 

QA  *  i  QQI  I2(  2  1-1)  *  RRI  I2(  2  1-1  )  I  /  2.000 
QB  ■  I  QQI  1212)1  *  RRII2I2H)  /  2.000 

Q  ■  (QA  ♦  QB  1/2.000 

ENO  IF 
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V-  1. 


. V.  ■ 


JUAb  1*. 


A 


ELSE  XF1CS0XR.ER.3 )  THEN 
C 

C  -  CONTACT  SURFACE  TRAVELING  LIFT  - 

C 

ZFIJ.ER.il  THEN 

RA  ■IQQA+RRA 1/2.000 
RB  ■  DQ*AA  ♦  RA 

CALL  SKJUNPI AA  ,RA ,SA , AR ,DQ ,VS ,G ,G1 ,H, AS ,QB ,$8 ) 

R  ■  RB 
GO  TO  200 

END  IF 

IF(X2(2I.ER.X2(1)1  THEN 

IF ( SIGMA ( 2  >1 )  .GE . SIGMA ( 1 >1 )  I  THEN 

<3  »  (QR( 1212)1  ♦  RR( 12(21)1  /  2.000 

ELSE 

R  ■  ( GQ(  Z2(  2  1-1 1  ♦  RR(  I2(  2  1-1 1 1  /  2.000 

END  IF 

ELSE  IF( SIGMA!  2,1  ).ER.  ( X2(  2  )-H ) )  THEN 
X2(  2  )  3  X2(  2  )  -  H 

12(2)  3  12(2)  -  1 

R  ■  ( QQ(  1 2 (  2  )  )«RR1  X2(  2  ) ) )  /  2.000 

ELSE 

QA  3  ( QQ(  1 2 (  2  )-l )  ♦  RR(  12(21-1)1  /  2.000 
RB  *  I  QQ(  12(2))  *  RR(  12(2)))  /  2.000 
R  «  (RA  ♦  RB  1/2.000 

END  IF 

END  IF 

200  SIGMA) 2,2)  3  OELT  *  R  ♦  SIGMA! 2,1) 

RCS*R 

C 

C  -  CALCULATE  EXPANSION  WAVE  SPEEO  - 

C 

TIME  «  TIME+OELT 
IF(J.ER.S)  THEN 

AW=(  <JQ(  ( N*1  )/2  )-RR(  ( N+l  )/2  ) )/(  2 .  DOO*S(  ( N*l  )/2 ) ) 

QH=Q 

IF( ( CSOIR  ) .EG. 2  )  THEN 

VHEW=-(  AW) 

VTEH*  QW-AW 

ELSE 

VHEWs  AW 
VTEH*  QW+AH 

END  IF 

END  IF 

IFIJ.EQ.3)  THEN 

SIGMA! 3.2  I  *  SIGMA! 3 >1)  ♦  VHEH*TIME 
SIGMA! 4 ,2  )  «  SIGMA! 4 >1 )  ♦  VTEW»TIME 
ELSE  IF! (CSOIR. ER. 2). ANO . I  SIGMA) 3 ,1 )  .GT . 0.000))  THEN 
SIGMA!  3.2)  «  SIGMA!  3.1)  ♦  VHEW*DELT 
SIGMA! 4 >2  )  »  SIGMAI 4,1  )  ♦  VTEW*0ELT 
ELSE  IF! ( CSOIR. EQ. 3  ) . ANO . I  SIGMA) 3,1  I . GT . 1 . 000  ) )  THEN 
SIGMA) 3 >2)  *  SIGMA! 3 .1)  ♦  VHEH«DELT 
SIGMA! 4,2  )  *  SIGMA! 4,1)  *  VTEH«DELT 

ENO  IF 
RETURN 
ENO 
C 

SUBROUTINE  SWEEP(N,H, SIGMA ,QQ,RR .S.DELT ,EE ,Q,A ,NEWRG,NENRR ,NEWS, 
CZ2,G2 , J.LNOOE , RNOOE .HALT ,LBNORY .LBOPRS , LBDPR , LBOTR , LBOOR ,QLBO ,G , 
CGI »RBNORY .RBOPRS ,RBOPR ,RBOTR .RBOOR ,RRBO ,BORY ,SK ) 

C 

C  MHHHHHHHHHHHHHHHHUHHHHHHHHHHHHHHHHHHMHUHHUHHHHW 


c  *  * 

C  »  SPACE  SWEEPING  SUBROUTINE  • 

C  *  * 

c 

c  -  VARIABLE  DEFINITIONS  - 

C 
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C  AAVC  -  AVERAGE  SPEED  OP  SCIJNO 

C  C NT ACT  -  S  DX6XT  VARIABLE  DENOTING  CONTACT  SURFACE 

C  LOCATION  ^DIRECTION  OF  TRAVEL*  ANO  IF  IT  CROSSES  A  NODE 

C  OELGQN  -  CHANGE  IN  G4  FROM  X  TO  X+l 

C  OELQGL  -  CHANGE  IN  GG  FROM  1-1  TO  X 

C  OELRRH  -  CHANGE  IN  RR  FROM  X  TO  X+l 

C  OELRRL  -  CHANGE  IN  RR  FROM  1-1  TO  X 

C  OELSH  -  CHANGE  IN  S  FROM  I  TO  X+l 

C  DELSL  -  CHANGE  IN  S  FROM  X-l  TO  X 

C  DELAH  -  CHANGE  IN  A  FROM  X  TO  I+l 

C  OELAL  -  CHANGE  IN  A  FROM  X-l  TO  X 

C  DELQH  -  CHANGE  IN  G  FROM  I  TP  1*1 

C  OELGL  -  CHANGE  IN  Q  FROM  1-1  TO  I 

C  DELX  -  INTERPOLATION  DISTANCE  ( LMD»OELT ) 

C  DLTA«h»  -  PREFIX  WHICH  INDICATES  THE  SPATIAL 

C  CHANGE  IN  **  FOR  ONE  TIME  STEP. 

C  XNTEGIK)-  RESULT  OF  INTEGRATING  ZIK) 

C  **INT  -  VALUE  OF  »»  INTERPOLATED  BETWEEN  NOOES 
C  ON  THE  CURRENT  TIME  LEVEL. 

C  LXX  -  NODE  DEFINING  THE  LEFT  INTERVAL 
C  «PRIM1K)-  SUFFIX  WHICH  INDICATES  THE  SPATIAL 
C  OERIVITIVE  OF  »  AT  THE  CURRENT  TIME  LEVEL. 

C  RXX  -  NODE  DEFINING  THE  RIGHT  INTERVAL 
C  SAVG  -  AVERAGE  ENTROPY 

C  SHOCK  -  S  DIGIT  VARIABLE  DENOTING  SHOCK 

C  LOCATION, DIRECTION  OF  TRAVEL,  AND  IF  IT  CROSSES  A  NODE 

C  wtSTEP  -  THE  CHANGE  IN  TIME  OF  M  AT  A  NODE 
C  USED  TO  STEP  UP  TO  THE  NEXT  TIME  LEVEL 

C  X  LOCATION  IN  SPACIAL  PLANE  1 1-1  >*H 

C  Z(K)  -  RIGHT  SIOE  OF  THE  K*TH  EQUATION. 

C 

INTEGER  I , RXX , LXX , SHOCK ,CNT ACT , 12 , J , LNODE , RNOOE .HALT , 

C  N , LBNORY , LBOPRS , RBNDRY .RBOPRS , SK , BDRY 

DIMENSION  SIGMA  I  4 , 2  )  ,S(  N I  ,Q(  N ) ,  A<  N )  ,121  4  I  ,QINT<  3  J  ,AINT(  3  )  ,Zf  3  ) , 
C  NEWQQt N I •NENRRI H ) ,NEWS( N I ,  INTEGI 3  )  .APRIMI 3 ) .QPRIMl 3  I , 

C  LNODE I 4  I .RNOOE ( 4 ) ,AAVG( 3  )  ,RR( N ) ,QQ( N  I 

DOUBLE  PRECISION  AAVG,SAVG,G2,X,H,SIGMA,QQ,RR,S,G,A,G,G1, 

C  OELQQH.OELQQL, OELRRH, OELRRL, DELSH, DELSL, 

C  OELAH.DELAL, DELQH, DELGL.DELT, 

C  QINT.AINT.QQINT, RRINT.SINT.EE, 

C  NEWQq .NEWRR .NEWS , LBDPR , LBDTR , LBDDR.QLBD , 

C  RRSTEP .SSTEP .INTEG.Z .OLTAQQ .QQSTEP , 

C  DLTARR .OLTAS .APRIM.QPRIM , 

C  RBOPR .RBOTR .RBDOR ,QRB0 

COMMON  AR.OQ.VS.W 
C 

C  -  COMPUTE  VELOCITY  ANO  SPEED  OF  SOUND  AT  EACH  NODE  - 

C 

DO  10  1=  l.N 

Q(  I  )  3  (QQ(I)  ♦  RR(  1 1  >  /  2.0000 

All)  3  (IQQIII  -  RRl  I ) )  /(  2.0000  *  SID)) 

10  CONTINUE 
C 

C  -  ADVANCE  LEFT  BOUNOARY  TO  NEW  TIME  STEP  - 

C 

BDRY  «  3 

IF!  121 1 ) .  EQ.  2  )  THEN 
SK  3  1 

ELSE  IF! SIGMA! 1,2  ).EQ. -2.D00 )  THEN 

SK  3  3 

ELSE 

SK  >  0 

END  IF 

CALL  BONORYI  Ql  1 ) ,Q(  2 )  >QLBD  ,A(  1 ) , A<  2  )  ,QQI  1 )  ,QQ(  2  )  ,RR1 1 )  ,RRI  2 ) , 

C  SI  1 ). SI  2), H, EE, OELT.LBNDRY, LBOPRS, LBDPR. LBDDR, LBDTR, 

C  J.NEWQQl 1 1.NEHRR1 1 1 .NEHSI 1 1 ,G,G1 ,G2 .HALT .BDRY ,SK  ) 

C 

C  -  AT  EACH  NODE  FROM  2  TO  N-l  DETERMINE  THE  BEST  ALGORITHM  TO  USE — 

C  -  TO  ADVANCE  THAT  NODE  TO  THE  NEXT  TIME  STEP  - 


157 


ooo  non  non  o  n  o 


C 

2-2 

11  XFI I.E6.N)  60  TO  1200 

X»FL0AT( 1-1  )»H 
DELQQH  ■  QQI 1*1 )  -  QQ(  I ) 

DELQQL  >  QQ(  I )  -  06(1-1 ) 

OELRRH  a  RR(I*1)  -  RR(  I ) 

OELRRL  ■  RRII)  -  RR(I-l) 

OELSH  a  SC  1*1)  -  SCI) 

DELSL  a  SCI)  -  SCI-1) 

DELAH  a  AC  1*1)  -  AC  I ) 

DELAL  a  All)  -  AII-1) 

OELQH  a  Q(I»1)  -  0(1) 

OELQL  -  0(1)  -  0(1-1) 

-  DEFINE  LEFT  SECTOR  AND  RIGHT  SECTOR  WRT  NOOE  EXAMINED— 

RXX  a  I  ♦  I 
LXX  ■  I 

-  TEST  FOR  SHOCK  - 

IF  (12(1). EQ. RXX)  THEN 
SHOCK  *  200 
GO  TO  20 

ELSE  IF  112(1). EQ. LXX)  THEN 
SHOCK  «  300 
GO  TO  20 

ELSE 

SHOCK  >  loo 

END  IF 
GO  TO  30 

-  DETERMINE  DIRECTION  SHOCK  IS  TRAVELING  - 

20  IF  (SIGMA(1,11.LT.SIGMA(1,2))  THEN 
SHOCK  «  SHOCK  ♦  20 

ELSE 

SHOCK  a  SHOCK  *  30 

END  IF 

-  DETERMINE  IF  SHOCK  CROSSES  A  NOOE  IN  THIS  TIME  INTERVAL 

IF  ( SHOCK. EQ. 220  )  THEN 

IF  ( SIGMA! 1 >2 ) . GE . ( X+H ) )  THEN 
SHOCK  a  SHOCK  ♦  l 

ELSE 

SHOCK  ■  SHOCK  ♦  2 

ENO  IF 

ELSE  IF  ( SHOCK. EQ. 230)  THEN 

IF  (SIGMAC 1,2). LE.X)  THEN 
SHOCK  a  SHOCK  ♦  1 

ELSE 

SHOCK  «  SHOCK  ♦  2 

END  IF 

ELSE  IF  (SHOCK. EG. 320)  THEN 

IF  ISIGMA(1,2).GE.X)  THEN 
SHOCK  a  SHOCK  ♦  l 

ELSE 

SHOCK  a  SHOCK  *  Z 

ENO  IF 

ELSE  IF  (SHOCK. EG. 330)  THEN 

IF  ( SIGMAC 1 > 2  ) . LE . ( X-H ) )  THEN 
SHOCK  ■  SHOCK  ♦  1 

ELSE 

SHOCK  a  SHOCK  ♦  2 

END  IF 

END  IF 
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c - test  for  contact  surface - 

c 

30  ZF  1X2(2). E6.RXX)  THEN 
CNTACT  ■  200 
60  TO  40 

ELSE  IF  1 121  2I.EQ.LXX)  THEN 
CNTACT  ■  300 
60  TO  40 

ELSE 

CNTACT  »  100 

END  ZF 
60  TO  50 
C 

C  - OETE RHINE  DIRECTION  CONTACT  SURFACE  IS  TRAVELING  - 

C 

40  IF  (SIGMA1 2 >1 ). LT .SIGMA)  2 >2 ) )  THEN 
CNTACT  ■  CNTACT  ♦  20 

ELSE 

CNTACT  ■  CNTACT  ♦  30 

END  IF 
C 

c  -  DETERMINE  IF  CONTACT  SURFACE  CROSSES  A  NODE  0URIN6  THIS  TIME  - 

C  -  INTERVAL  - 

C 

IF  (CNTACT.EQ.220)  THEN 

IF  (SIGMAI 2>2  )  .GE . I X+H ) )  THEN 
CNTACT  ■  CNTACT  ♦  1 

ELSE 

CNTACT  ■  CNTACT  ♦  2 

END  IF 

ELSE  IF  ( CNT ACT . EG . 230 )  THEN 

IF  ( SIGMAI 2  >2  ) . LE .X )  THEN 
CNTACT  a  CNTACT  ♦  1 

ELSE 

CNTACT  *  CNTACT  ♦  2 

END  IF 

ELSE  IF  ( CNT ACT. E6. 320)  THEN 

IF  l SIGMAI 2  >  2  ) .  GE .  X )  THEN 
CNTACT  »  CNTACT  ♦  1 

ELSE 

CNTACT  «  CNTACT  ♦  2 

END  IF 

ELSE  IF  ( CNTACT. Eq. 330)  THEN 

IF  (  SIGMAI  2  >  2  ) . LE . ( X-H  )  )  THEN 
CNTACT  a  CNTACT  ♦  1 

ELSE 

CNTACT  ■  CNTACT  ♦  2 

END  IF 

END  IF 
C 

C  - CHECK  IF  EITHER  A  SHOCK  OR  CONTACT  SURFACE  HITHIN  H  OF  THIS  NODE- 

C  - OETERMINE  PROPER  ALGORITHM  TO  USE  FOR  CALCULATING  EXTENDED  REIMAN- 

C  - VARIABLE  CHANGE  ALONG  CHARACTERISTICS  AT  THIS  NODE - 

C 

50  IF  ( SHOCK. EQ. 100. OR. CNT ACT. Eq. 100)  THEN 
C 

C  - NEITHER  A  SHOCK  NOR  A  CONTACT  SURFACE  EXIST  NEAR  NODE - 

C 

IF  (SHOCK. Eq. 100. ANO. CNTACT. Eq. 100)  THEN 
C 

C  - TEST  FOR  SUBSONIC  OR  SUPERSONIC  FLOW--- 

C 

IF  (DABS! q<  I )  I.LT.AI  I ) )  THEN 

IF(  I  .Eq.  ( I2(  1  )-2  ) )  THEN 
XFM.Eq.H  THEN 

IFl SIGMA) 1,1  ).LE. SIGMA) 1,2)) 
GO  TO  200 

END  ZF 

END  IF 


THEN 
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END  IF 

XF(X.EQ.(X2(1)+1))  THEN 
IF(J.EQ.l)  THEN 

IF( SIGMA 1 1 ,1 ) .GE .SIGMA( 1 >2  ) )  THEN 
GO  TO  300 

END  IF 

END  IF 

END  IF 
GO  TO  100 

ELSE 

IF  ( SIGMA! 2.1). LE .SIGMA! 2,2))  THEN 
GO  TO  200 

ELSE 

GO  TO  300 

END  IF 

END  IF 

EM)  IF 

- SHOCK  OR  CONTACT  SURFACE  ON  LEFT,  HEADED  RIGHT,  NO  NODES  CROSSED — 

IF  ( SHOCK . EG . 322 . OR . CNT ACT . EQ . 32 2 )  THEN 
IF  (OABSIQII n.LT.AII))  THEN 
GO  TO  300 

ELSE 

GO  TO  500 

END  IF 

END  IF 

- SHOCK  OR  CONTACT  SURFACE  ON  LEFT,  HEADED  LEFT,  NO  NODES  CROSSED — 

IF  ( SHOCK. EQ. 332. OR. CNT ACT. EQ. 332)  THEN 
GO  TO  300 

END  IF 

- SHOCK  OR  CONTACT  SURFACE  ON  LEFT,  HEADED  RIGHT,  NODE  IS  CROSSED— 

IF  (SHOCK. EQ. 321. OR. CNTACT.EQ. 321)  THEN 
GO  TO  400 

ENO  IF 

- SHOCK  OR  CONTACT  SURFACE  ON  LEFT,  HEAOED  LEFT,  NODE  IS  CROSSED - 

IF  (SHOCK. EQ. 331. OR. CNTACT.EQ. 331)  THEN 
GO  TO  300 

ENO  IF 

- SHOCK  OR  CONTACT  SURFACE  ON  RIGHT,  HEAOED  RIGHT,  NO  NODES  CROSSED — 

IF  (SHOCK. EQ. 222. OR. CNTACT.EQ. 222)  THEN 
GO  TO  200 

ENO  IF 

- SHOCK  OR  CONTACT  SURFACE  ON  RIGHT,  HEADED  LEFT,  NO  NODES  CROSSED— 

IF  (SHOCK. EQ. 232. OR. CNTACT.EQ. 232)  THEN 
IF  (DABSIQ(II).LT.A(I))  THEN 
GO  TO  200 

ELSE 

GO  TO  500 

END  IF 

ENO  IF 

—SHOCK  OR  CONTACT  SURFACE  ON  RIGHT,  HEADED  RIGHT,  NODE  CROSSED - 

IF  (SHOCK. EQ. 221. OR. CNTACT.EQ. 221)  THEN 
GO  TO  200 

ENO  IF 

- SHOCK  OR  CONTACT  SURFACE  ON  RIGHT,  HEADED  LEFT,  JUMPS  NODE - 
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GO  TO  400 
ENO  IF 

- BRANCH  HERE  IF  SHOCK  TO  RIGHT  OF  CONTACT  SURFACE - 

- DETERMINE  PROPER  ALGORITHM  TO  USE  FOR  CALCULATING  EXTENDED  REIMAN— 

- VARIABLE  CHANGE  ALONG  CHARACTERISTICS  AT  THIS  NOOE - 

IF  (SIGMA( 1»1 ).GT.SIGMA( 2 >1 ) )  THEN 

- SHOCK  ON  LEFT  ,  HEADED  RIGHT ,  NO  NOOE  JUMPED - 

IF  (SHOCK. EG. 322)  THEN 

IF  ( CNTACT.EQ.322  )  THEN 

IF  (  OABSI  Q( I )  ) . LT . A( I ) )  THEN 
GO  TO  300 

ELSE 

GO  TO  500 

END  IF 

ELSE  IF  (CNTACT.EQ.332)  THEN 
GO  TO  800 

ELSE  IF  ICNTACT.EQ.331)  THEN 
GO  TO  800 

ELSE 

GO  TO  900 

EM)  IF 

END  IF 

- SHOCK  ON  LEFT,  HEAOEO  LEFT,  NO  NOOE  JUMPED - 

IF  (SHOCK. EG. 332)  THEN 

IF  (CNTACT.EQ.322)  THEN 
GO  TO  300 

ELSE  IF  (CNTACT.EQ.332)  THEN 
GO  TO  600 

ELSE  IF  (CNTACT.EQ.331 )  THEN 
GO  TO  600 

ELSE  IF  (CNTACT.EQ.321)  THEN 
GO  TO  700 

ELSE 

GO  TO  900 

END  IF 

END  IF 

- SHOCK  ON  LEFT,  HEAOEO  RIGHT,  JUMPS  NODE - 

IF  (StKXK.EQ.321 )  THEN 

IF  (CNTACT.EQ.322.OR.CNTACT.EQ.321)  THEN 
GO  TO  400 

ELSE  IF  ( CNTACT.EQ.332. OR. CNTACT.EQ.331)  THEN 
GO  TO  800 

ELSE 

GO  TO  900 

END  IF 

END  IF 

- SHOCK  ON  LEFT,  HEADED  LEFT,  JUMPS  NOOE - 

IF  ( SHOCK. EQ. 331)  THEN 

IF  ( CNTACT.EQ.322. OR. CNT ACT. EQ. 321 )  THEN 
GO  TO  700 

ELSE  IF  (CNTACT.EQ.332. OR. CNTACT.EQ.331)  THEN 
GO  TO  600 

ELSE 

GO  TO  <*00 

END  IF 

END  IF 

- SHOCK  ON  RIGHT,  HEAOEO  RIGHT,  NO  NOOE  JUMPED— 
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IF  ( SHOCK. ES. 222)  THEN 

IF  (CNTACT.EQ.222)  THEN 
SO  TO  200 

ELSE  IF  (CNTACT.ES. 232. OR. CNTACT.ES. 251)  THEN 
GO  TO  800 

ELSE  IF  ICNTACT. EQ.321.0R.CNTACT.ES. 522)  THEN 
GO  TO  400 

ELSE  IF  (CNTACT.ES.332.OR.CNTACT.EQ. 351)  THEN 
GO  TO  800 

ELSE 

GO  TO  900 

END  IF 

END  IF 
C 

C - SHOCK  ON  RIGHT,  HEADED  LEFT,  NO  NODE  JUMPED— 

C 

IF  ( SHOCK. EG. 232)  THEN 

IF  (CNTACT.EQ.222)  THEN 

IF  ( SIGMA ( 1 ,2  I . LT . SIGMA) 2,2))  THEN 
GO  TO  700 

ELSE 

GO  TO  200 

END  IF 

ELSE  IF  ( CNTACT.EQ.231.OR.CNTACT.EQ. 232)  THEN 
GO  TO  600 

ELSE  IF  (CNTACT.ES. 221)  THEN 
GO  TO  700 

ELSE  IF  (CNTACT.ES. 322)  THEN 
GO  TO  400 

ELSE  IF  (CNTACT.ES. 332. OR. CNTACT.ES. 331)  THEN 
GO  TO  600 

ELSE  IF  ( SIGMA! 1,2). LT. SIGMA ( 2,2))  THEN 
GO  TO  710 

ELSE 

GO  TO  400 

END  IF 

END  IF 
C 

C  - SHOCK  ON  RIGHT,  HEADED  RIGHT,  JUMPS  NODE--- 

C 

IF  ( SHOCK.ES- 221 )  THEN 

IF  ICNTACT. EQ. 221. OR. CNTACT.EQ. 222 )  THEN 
GO  TO  200 

ELSE  IF  ICNTACT. EQ. 251. OR. CNTACT.ES. 252 >  THEN 
GO  TO  300 

ELSE  IF  ICNTACT. ES.321.OR.CNTACT.ES. 322)  THEN 
GO  TO  400 

ELSE 

GO  TO  800 

ENO  IF 

END  IF 
C 

C  —SHOCK  ON  RIGHT,  HEADED  LEFT,  JUMPS  NOOE  — 

C 

IF  ICNTACT. EQ.221.0R.CNTACT.es. 222)  THEN 
GO  TO  700 

ELSE  IF  (CNTACT.EQ. 231. OR. CNTACT.EQ. 232)  THEN 
GO  TO  600 

ELSE  IF  (CNTACT.ES. 322)  THEN 

IF  (SIGMA! 1,2  I.LT.SIGMAI 2,2 ) )  THEN 
GO  TO  710 

ELSE 

GO  TO  400 


ELSE 

60  TO  800 

ENO  IF 
C 

C  - BRANCH  HERE  IF  SHOCK  IS  TO  LEFT  OF  CONTACT  SURFACE - 

C  - DETERMINE  PROPER  ALGORITHM  TO  USE  FOR  CALCULATING  EXTENDED  REIMAN- 

C  - VARIABLE  CHANGE  ALONG  CHARACTERISTICS  AT  THIS  NOOE - 

C 

ELSE  IF  ( SIGMA! 1 > 1 ) . LT . SIGMAI  2,1))  THEN 
C 

C - SHOCK  ON  LEFT,  HEADED  RIGHT,  NO  NOOE  CROSSED— 

C 

IF  (SHOCK. EO. 522  I  THEN 

IFtCNTACT.EQ. 322 . OR . CNTACT . EQ. 121  )  THEN 
GO  TO  oOO 

ELSE  IF  I CNTACT.EQ. S3 2)  THEN 

IF  ( SIGMAI 1 ,2  I ■ GT . SIGMAI 2,2)1  THEN 
GO  TO  700 

ELSE 

GO  TO  300 

ENO  IF 

ELSE  IF  I CNTACT. EG. 331 1  THEN 
GO  TO  700 

ELSE  IF  (CNTACT. EG. 222. OR. CNTACT. EO. 221 )  THEN 
GO  TO  600 

ELSE  IF  I CNTACT. EG. 232)  THEN 
GO  TO  400 

ELSE  IF  (SIGMAI  1,2  J.GT. SIGMAI  2,2))  THEN 
GO  TO  710 

ELSE 

GO  TO  400 

ENO  IF 

ENO  IF 
C 

C  —SHOCK  ON  LEFT,  HEADED  LEFT,  NO  NOOE  CROSSED— 

C 

IF  I  SHOCK. EO. 332)  THEN 

IF  ( CNTACT . EQ. 322 )  THEN 
GO  TO  800 

ELSE  IF  ( CNTACT . EQ. 332  )  THEN 
GO  TO  300 

ELSE  IF  I CNTACT . EQ. 321 )  THEN 

GO  TO  aoo 

ELSE  IF  (CNTACT. EQ. 222. OR. CNTACT. EQ. 221  I  THEN 
GO  TO  300 

ELSE  IF  (CNTACT. EO. 231. OR. CNTACT. EQ. 232  I  THEN 
GO  TO  400 

ELSE 

GO  TO  400 

ENO  IF 

ENO  IF 
C 

C  —SHOCK  ON  LEFT,  HEADED  RIGHT,  JUMPS  NOOE  — 

C 

IF  ( SHOCK. EQ. 321)  THEN 

IF  (CNTACT.EQ. 322. OR. CNTACT. EQ. 321 1  THEN 
GO  TO  600 

ELSE  IF  ( CNTACT.EQ. 332. OR. TNT ACT. EQ. 331  I  THEN 
GO  TO  710 

ELSE  IF  (CNTACT. EQ. 222. OR. CNTACT. IQ. 221  )  THEN 
GO  TO  600 

ELSE  IF  ( CNTACT.EQ. 231 )  THEN 
GO  TO  710 

ELSE  IF  (SIGMAI 1,2). GT. SIGMAI 2,2  I)  THEN 
GO  TO  710 

ELSE 

GO  TO  400 

ENO  IF 

ENO  IF 
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-shocx  on  left ,  hcamd  left.  jumps  non 


IP  I  SHOCK  . EO. 331)  THCN 

IP  (CNTACT . EG.  322 .OR .CMT ACT . EO. 321 )  THCN 
00  TO  MO 

f LSI  IP  I CNTACT. EG. 331. OR. CUT ACT. EG. 332  I  THEN 
CO  TO  100 

I  LSI  IP  (CNTACT. EG. 221. OR. CNTACT. EG. 222)  THCN 
00  TO  000 

I  LSI 

00  TO  400 

ENO  IP 

I  NO  IP 

SHOCK  ON  RIGHT ,  HI  AMO  RIGHT,  NO  NON  CROSSED— 

IP  I  SHOCK  . 10.222  I  THEN 

IP  ICNTACT.E0.222.0R.CNTACT.E0  221  I  THEN 
GO  TO  *00 

ELSE  IP  (CNTACT  EG. 231  I  THEN 
GO  TO  UO 

ELSE  IP  CONTACT  EG. 212  I  THEN 

IP  I  SIGMA! 1 ,2  1 . 6T . SIGMA I  2,211  THEN 
GO  TO  700 

ELSE 

GO  TO  200 

END  IP 

ELSE 

00  TO  *00 

ENO  IP 

ENO  IP 

-SHOCK  ON  RIGHT.  HEAMO  LEFT,  NO  NOMS  CROSSED— 

IP  (SHOCK  EO.  232  I  THEN 

IP  CCNTACT  EO  222. OR  CNTACT  EO  221  I  THEN 
GO  TO  tOO 

ELSE  IP  i CNTACT  EO. 232  l  THEN 

IP  (OAGSiQC  I  M  LT  .  Ai  I  M  THEN 
GO  TO  200 

ELSE 

00  TO  soo 

ENO  IP 

ELSE 

GO  TO  TOO 

EM)  IP 

ENO  IP 

-SHOCK  ON  RIGHT.  HEADED  RIGHT.  JUMPS  NOOE — 

IP  (SHOCK  EO. 221  I  THCN 

IP  (CNTACT  EO  221. OR. CNTACT  EO  222 »  THEN 
GO  TO  »00 

ELSE  IP  (CNTACT  EO. 231 (  THEN 
GO  TO  710 

ELSE  IP  *  CNTACT  EO  232  I  THEN 
80  TO  700 

ELSE 

00  TO  TOO 

ENO  IP 

ENO  IP 

•SHOCK  ON  RIGHT.  HCAOED  LEFT,  JUMPS  NODE  — 

IP  CNTAC'  £9  :zi  JR  CNTAC'  EO . ZZZ  'HEM 
GO  TO  ROO 

ELSE  IP  (CNTACT  to. 231. OR. CNTACT  EO  232 >  THEN 
GO  TO  400 


I 


UOOU  u  u  u  u 


else 

00  TO  <00 

cmo  if 

- SHOCK  AND  CONTACT  SUMAC!  AM  AT  TNI  SAMI  LOCATION  AFTIR  TIMS - 

— zero — 


■  LSI  iril SHOCK. 19. Ztt I. AND. <CNT ACT. EQ. 222 DTHCN 
SHOCK  ■  S21 
CNTACT  •  S2X 
00  TO  <00 

ELSE  IFI <  SHOCK . EQ. J22  > . ANO. I CNTACT . EQ. J22  1 1THEN 
IKIOASSIOII  ll.LT.AII  I)  THEN 
CO  TO  SOO 

ELSE 

60  TO  SOO 

ENO  If 

ELSE  IPl (SHOCK. E9.SJ2I. ANO. (CNTACT. f0.JS2»>  THEN 
SHOCK  •  2J1 
CNTACT  ■  211 
CO  TO  <00 

ELSE  Ifl  (SHOCK.  EQ.232I.  ANO.  (CNTACT.  EQ.2S21I  THEN 
IF(0A93(Q(  X  I  I.LT.AII  I)  THEN 
CO  TO  200 

ELSE 

90  TO  SOO 

eno  xr 

ELSE 

GO  TO  720 

ENO  If 

- CALL  CQNOITION  SUOMUTIM  MHICH  CONTAINS  ALGORITHM  THAT  IS— 

— NUHERICALLY  BEST  SUITED  to*  THE  SITUATION  AT  NOOE  I  — 


C 


c 


c 


c 


c 


c 


100  CALL  CONOKOI  I  l.9<  1*1  l.AI  I  >.A(I*l  t.RRi  I  l.M(  I  l.SI  I  I.DELML, 

C  Of  LRRH  .OELSH  .OELSL  .OELQM  .OELAH  .0EL9L  .OELAL  .H  .EE  . 

C  0tL7.  OOINr.ARINT.SXNT.OPAIN,  AMIN,  MlffT. 91  I -1  ),AJ  I-])) 

00  TO  1000 

200  CALL  C0N02(  Q(  I  1 .9*  I  *1  I  >A(  I  I  >A(  1-1  I  .RSI  I  >  .QQI  I  I  ,S(  I  1 .0EL09L  > 

C  OELRRL  .OELSL  .Of  L9L  .OELAL  .OELT  ,H,EE  .MINT  , RHIWT  ,SINT  , 

C  '3PRIN.APRIN.AINT  I 

00  TO  1000 

SOO  CALL  CONOSiOl  I  1,91  1*1  >.AI  I  ),A(I*1  I.RRl  I  I.Ml  I  I.Sl  I  I.OELQQH, 

C  Of  l RRH  .OELSH  , 0EL9H.DE L AH  .OELT  ,H.EE  .MINT  .RRINT  .SINT  , 

C  9PRIH.AP0IH.AINT I 

00  TO  1000 

<00  CALL  CONO<< I .SHOCK .CNTACT .J.LNOOE .RNOCE .99STEP. RRSTEP. SSTEP ( 

00  TO  1100 

SOO  CALL  CaN0f(0<!l.9«I'll.9(I«ll.AIX>,  AII-11,  Al  I  *1  >  ,RR(  I  >  ,M(  I  > . 

C  SI  I  I.OELML  .Of  LOOM. OELRRL  .Of  LRRH ,0f  LSI  , OELSH, Of  LOL  . 

C  Of  LOH, OELAL  ,0f  LAH.H.ff  .OELT  .MINT  .RRINT  ,SINT  , AINT  , 

C  9PR1H. APRIH. SHOCK .CNTACT  I 

oo  to  looa 

<00  CALL  C»«*I  SHOCK,  CNTACT,  HALT  I 
OOSTf P*0  00 
RRST|»«0  00 
SSTEP*0  00 
00  TO  UOO 


’00  CALL  CONO?«  SHOCK  .CNTAC’ .OELT  .3ISNA  .12  ,  1  ,H  .malt  ,01  I  I 
p’  i  sn 

RRST|P*0  00 
SSTf P»0  00 
90  TO  1100 


165 


£ 
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i 

*g 


i 
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CALL  C0ND7NI SMOCK  .CNTACT ,01 LT .SIGMA ,1* ,1 ,H .HALT ,«t I )  ,X ) 
99STEP*0.00 
RRSTEP*0.00 
SSTEMO.DO 
GO  TO  UOO 


CALL  C0ND7SI SIGMA, HALT, SHOCK, CNTACT) 
GGSTEP«0.00 
BASTE P *0.00 
SSTEP-0.00 
GO  TO  1100 


CALL  CONOBI SHOCK .CNTACT, HALT) 
99STEP»0 . 00 
RRSTEP*0.00 
SSTEP*0 . 00 
GO  TO  1100 


PRINT  •  , 'AN  IMPOSSIBLE  SITUATION  EXISTS  -  ERROR* 
QQSTEP-O.DO 
RRSTEP*0. 00 
SSTEP«0.00 
HALT  •  1 
GOTO  1100 


CALCULATE  OLTA  99,  OLTA  RR  A  OLTA  S 


OLT A99*Q9IMT -9Q( I ) 

OLTARR*RRINT-RRI I ) 
0LTAS>S1NT-SI I  I 


CALCULATE  ZIKI'S 


AAVGI  1  )a(  AINTI 1  )«A<  I )  1/2.0000 
AAV9I  3  )■»  AINTI  3  >*A<  I  )  1/2.0000 
AAVGI 2  1*0.0000 
SAVQ  a  ( SINT *S(  I  I  1/2 . 000 

Zll  >•-!  1 . 0000/G2  )«AAVGI  1  )»f  SAVG-G2  M»l  9PRIMI 1  )-GZ*APRIHI  1 ) ) 
Z(  3  )»l  1 . 0000/02  l*AAVG(  3  )•(  SAVG-G2  )»l  9PRIMI  3  )*GZ*APRIMI  3  )  ) 
Zl  Z  1*0 . 0000 


INTEGRATE  THE  Z(K)'S 


INTEGI Z  1*0 . 0000 
INTEGI  1  )*ZI  1  l*0ELT 
INTEGI  3  )*ZI  3  l*DELT 


SOLVE  THE  E9UATX0N 


99STEP*0LTAQQ+ INTEGI  1 ) 
RRSTEP*OLTARR* INTEGI  3  I 
SSTEP*OLT AS* INTEGI  2  I 


STORE  THE  SOLUTION 


NEMQ9I I  1*001  I  l*QQSTEP 
NEMRRI I  l*RRI  I  HRRSTEP 
NEMSI I  l*S<  I  l*SSTEP 


SO  TO  NEXT  NODE 


I»I*1 
GOTO  11 
CONTINUE 

•ORY  ■  2 

IF!  121  1  I.E9  N)  THEN 
SK  «  1 


•;**i 


t£SA 


ELSE  IF (SIGMA! 1 ,2 ) .  E9.3 .000  I  THEN 
SK  «  3 

ELSE 

SK  •  0 

ENO  IF 

CALL  BONORYI 61  N  )  ,61  N-l )  >QRBO  »AI  N )  »AI  N-l )  ,96<  N )  ,661 N-l )  ,RRI  N  )  , 
C  RRI N-l > ,SI  N ) ,SI N-l  I  .H ,EE ,DELT .RBNDRY , RBDPRS ,RBOPR ,RBDDR , 

C  RBOTR , J ,NEH6QI N )  .NEHRRI N ) ,NEHS( N  > ,6,61 ,62 ,HALT ,BORY ,SK ) 


UPDATE  THE  VARIABLES 


I-l 

IF  ( I . EQ.N+1 )  GOTO  1220 
RRID-NEHRRII) 

QQI I  )»NEWQGI I ) 

S(  I  )-NEWS(  I  ) 

1-1*1 
GOTO  1210 
CONTINUE 
RETURN 
ENO 


CONOITION  SUBROUTINES  1-B 


SUBROUTINE  CONOll 61 ,«IP1 ,AI , AIP1 ,RR ,60 .S .0EL96L .OELRRH .DELSH . 

C  OELSL  ,OELQH,OELAH,OELQL .OELAL .H .EE .OELT .60INT . 

C  RRINT  .SINT  ,QPRIM,APRIM,AINT  .4IM1.AIM1  I 


SUBROUTINE  CONOITION  1 


USED  WHEN  NO  CONTACT  SURFACES  NOR  SHOCKS  MITHIN  H  OF  NODE - 

AND  FLOW  IS  SUBSONIC - 

CALCULATES  OQINT .RRINT .SINT .qPRIM.APRIM - 

USES  FORHARO  ANO  BACKWARD  DIFFERENCE  SCHEMES - 


VARIABLE  DEFINITIONS 


AI  -  All) 

AIM1  -  AII-1)' 

AIP1  _  AI 1*1 ) 

El K I  -  ACTUAL  ERROR  IN  CHARACTERISTIC  SLOPE  CALCULATION. 
LMO  -  SLOPE  OF  THE  CHARACTERISTICS  (0*A,8,6-A) 

61  -  611) 

9IM1  -  911-1) 

6IP1  -  611*1) 


DIMENSION  LMOI  3 ) ,OELXI  3 ) ,9INT(  3 ) ,AINTI 3 ) ,E(  3 ) .6PRIMI  3 ) .APRIMI  3  I 
INTEGER  K 

DOUBLE  PRECISION  91 ,OIPl .AI ,AIP1 ,RR .66 .S . AIM1 ,GIM1 , 

C  DEL66L .OELRRH .DELSH .OELSL .0EL6H .0ELAH.DEL6L . 

C  OELAL .H.EE .OELT ,LMO ,DELX .9INT , AINT ,E , 

C  GQINT, RRINT, SINT, 6PRIM.APRIM 


-INITIAL  ESTIMATE  OF  CHARACTERISTIC  SLOPES- 


LM0I1)  ■  6IM1  *  AIM1 

LMOI  2  )  *  61 

LMD(J)  *  3IP1  -  AIP1 


- -CALCULATE  LINEARLY  INTERPOLATED  VALUES  OF  6  ANO  A- 


ononnocioon  o  non  onnn  non  n  noon 


c 

10  It  ■  1 

20  If  IK. IT. 4  )  THEN 

DC  LX  I K  I  •  DCLT  »  LMDIK) 


IP  ILMDIM.LT. 0 

.00001 

THEN 

QINTI  M 

a 

QI  - 

1 OCLXI M 

» 

OELQH 

/ 

H) 

AINTI  K  1 

a 

AI  - 

1 OCLXI M 

» 

DELAH 

/ 

HI 

ELSE 

QINTI  M 

a 

QI  - 

IDELXIK) 

• 

DELQL 

/ 

H) 

AINTI  M 

a 

AI  - 

IDELXIM 

« 

OELAL 

/ 

H  1 

END  IP 
K  ■  K  ♦  1 
GO  TO  20 

END  IP 

- -CALCULATE  ERROR  BETWEEN  ESTIMATED  SLOPE  AND  NEW  SLOPE - 

- PROM  MEN  INTERPOLATED  VALUES - 

till  ■  OABSILMOI1I  -  IQINTUI  «  AIMTIim 
El  2)  ■  0ABSILM0I2)  -  IQXNT(2)I) 

CIS)  >  DABSILMOIS)  -  I QINTI 3  I  -  AINTI III) 

LMDIl)  ■  QINTI  1 )  «  AINTI  1 1 

LMDI2)  •  QINTI  2  I 

LM0I3)  ■  QINTI  3)  -  AINTI  5 ) 

- -COMPARE  ERROR  TO  ERROR  TOLERANCE  LEVEL.  ITERATE  TIL  MET - 

IP  IIEI1I.GT.EEI.OR.IE<2>.GT.EE).OR.1EIS).GT.EE)>  00  TO  10 

- CALCULATE  LINEARLY  INTERPOLATED  VALUES  OP  REIMAN - 

- VARIABLES  ANO  NOOIF1CO  ENTROPY  AT  POINT  A - 

OQINT  ■  QQ  -  I0CLXI1I  •  i OELQQL  /Hi) 

IP  I LMDI 2  I . LC . 0 . 0000  I  THEN 

SINT  •  S  -  I OCLXI 2  I  ■  I OELSH  /  HI) 

ELSE 

SINT  «S-  I OCLXI 2  I  •  I0CLSL  /  HI) 

ENO  IP 

RRINT  ■  RR  -  I OCLXI 3 1  ■  I0CLRRH  /  Hil 
—  CALCULATE  SPATIAL  DERIVATIVES  — 


QPRIMI  1  I 
QPRIMI 2  I 
QPRIMI 3  I 
APRIMI 1  I 
APRIMI 2  I 
APRIMI 3  I 
RETURN 
ENO 


l  QI-QIM1  l/H 
0.D00 

<  QIPl-QI  l/H 
I  AI-AIM1  l/H 
0.000 

I  AIPI-AI  l/H 


SUBROUTINE 

C 

C 


CON02I QI .QIM1 i AI , AIM1 .RR ,QQ ,S, OELQQL .DCLRRL .OELSL . 
OCLQL .OCLAL  > DC LT .H.EE .QQINT .RRINT ,SINT .QPRIM. 
APRIM.AINT I 


SUBROUTINE  C0N0ITI0N  2 


■ 


—  BACKWARD  DIFFERENCE  ALGORITHM  FOR  CALCULATING  — 

--  RRINT .QQINT, SINT .APRIM.QPRIM.AINT  --- 

DIMENSION  LMOi  1  I.OEIXI  I  I.GINTl  3  I, AINTI  3  I .  E  ■  3  l. QPRIMI  3  l. APRIMI  I  I 
INTEGER  K 

DOUBLE  PRECISION  QI .QIM1 ,AI ,AIM1 ,RR ,QQ,S .OELQQL .DCLRRL .OELSL . 

C  OCLQL,  OELAL.DELT.H,  EE  .AINT  ,QINT  .LMO.DELX.l , 
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C  POINT .RRINT, SINT. OPRIM, APRIM 

C 

C  - INITIAL  ESTIMATE  Of  CHARACTERISTIC  SLOPES - 

C 

LMOI 1 )  a  0IM1  ♦  AIH1 
LMOI  S)  ■  OX 
LMOI  SI  ■  01  -  AX 

c 

C  - -CALCULATE  LINEARLY  INTERPOLATED  VALUES  Of  0  AND  A - 

C 

10  K  ■  1 

20  Xf  (K.LT.A)  THEN 

OELXIK)  ■  DELT  *  LMOIK) 

9XNTIK)  ■  01  -  I DC  LX  I K  I  *  OELQL  /  H) 

AINTI  K  I  ■  AX  -  ( OELXI  K  )  R  DELAL  /  HI 
K  «  K  «  1 
60  TO  20 

END  Xf 
C 

c  - -CALCULATE  ERROR  BETNEEN  ESTIMATED  SLOPE  AND  NEH  SLOPE-  --- 

C  - FROM  NEH  INTERPOLATED  VALUES - 

C 

Ell)  •  OABS(LMOIl)  -  IQINTIl)  ♦  AINTI 1)1) 

Bill  ■  DABS  I  LMOI  2  )  -  OXNTI  2  1 1 

SIS)  a  OABSILMDISI  -  IOINTISI  -  AINTI  311) 

C 

LMOI  II  •  OXNTI  1 1  ♦  AINTI  1) 

LMDI2I  a  QINTIXI 
LMOI SI  «  QINTIS)  -  AINTI SI 
C 

c  - COMPARE  ERROR  TO  ERROR  TOLERANCE  LEVEL.  ITERATE  TIL  MET - 

C 

If  IIEIll.6T.EEI.0R.lEt2l.eT.EEI.0R.lEISI.6T.EEn  00  TO  10 
C 

c  - CALCULATE  LINEARLY  INTERPOLATED  VALUES  Of  REIMAN - 

C  - VARIABLES  ANO  MODIFIED  ENTROPY  AT  POINT  A - 

C 

POINT  •  00  -  I  DELHI  1  I  a  tOELQOL  /  Hit 
SINT  as-  I  DELHI  2  I  a  I  OELSL  /  Nil 
RRINT  a  RR  -  I  DELHI  SI  a  I OELRRL  /  HI) 

C 

C  —  CALCULATE  SPATIAL  DERIVATIVES  — 

C 

OPRIMI  1  I  a  I  QI-6IM1  l/H 
OPRIHI 2  I  >  0 . 000 
OPRIMI  SI  «  (  QX-OIMX  l/H 


APRIMI  1 1  ■  IAI-AIM1I/H 
APR1MI 2  I  ■  0.000 
APRIMI  SI  a  I  AI-AIM1  l/H 
RETURN 
ENO 

SUBROUTINE  CONOSI 01 .0XP1 .AX ,AXP1 ,RR ,00 ,S .DE LOOM .DELRRH .DELSH , 
C  OELOHiDE LAM, DELT ,H.EE .OQINT , RRINT .SINT .QPRIM, 

C  APRIM, AINTI 


SUBROUTINE  COrCITXON  S 


fOPMARO  DIFFERENCE  ALGORITHM  FOR  CALCULATIN6  — 
RRINT , POINT, SINT, APRIM, 6PRIM, AINT  --- 


DIMENSION  .MCI  1  i. OELXI  1  l.OINTI  1  i, AINTI  5  I. El  5  I. OPRIMI  I  I, APRIMI  J  I 
INTEGER  * 

DOUBLE  PRECISION  01 .QIP1 ,AI ,AIP1 ,RR ,00 ,S , 

C  DELOOM, DELRRH, Of LSM, DELON, DELAH, 
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c 

c 


DCLT  ,H  ,CC  .AXNT  ,OXNT  .LMO.DELX.E, 
MINT  (MXNT  .SINT  ,4WN  .APRIM 


C 

C 

c 


c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


- INITIAL  ESTIMATE  Of  CHARACTERISTIC  SLOPES - 

LMO(l)  ■  OX  ♦  AX 
LMDIE)  ■  OX 
LMD(J)  a  OXP1  -  AXP1 

- -CALCULATE  LINEARLY  INTERPOLATED  VALUES  Of  0  AND  A- 

10  K  a  x 

20  If  IK.LT.4)  THEN 

OELXIKI  *  CELT  a  LMOI K  I 

QINTIKI  ■  01  -  I OELXIK )  •  0EL9H  /  HI 

AINTIK)  •  AI  -  ( OELXIK)  •  OELAH  /  H) 

K  a  K  ♦  I 
00  TO  20 

END  If 


- -CALCULATE  ERROR  SETNEEN  ESTIMATED  SLOPE  AND  NEM  SLOPE - 

- PROM  NEM  INTERPOLATED  VALUES - 

1(11  •  OASSILMOIXI  -  I0INTI1I  *  AINT(  111) 

Bill  a  0ASSILM0I2)  -  0INT<2>) 

EISI  ■  DASSILMOm  -  (OINTIS)  -  AINT(JU) 

LMOI  1 )  ■  0INTI1)  ♦  AINTI  1 1 

LMOI  2  I  «  QINTI2I 

LMOI SI  ■  OINTIS  I  -  AINTI SI 


- COWARE  ERROR  TO  ERROR  TOLERANCE  LEVEL.  ITERATE  TIL  MET - 

XfllElll.CT.EEI. ON. IEI2). 0T.EE).OR.IEISI.6T.EE>>  60  TO  10 

- CALCULATE  LINEARLY  INTERPOLATED  VALUES  OP  REIMAN - 

- VARIASLES  AND  MODIFIED  ENTROPY  AT  POINT  A - 

MINT  a  M  -  I DELXI 1 1  a  DELOOM  /  Hi 
SINT  as-  I  OELXI  2  I  a  DELSM  /  HI 
PRINT  ■  RR  -  I OELXI SI  •  OELRRH  /  HI 

—  CALCULATE  SPATIAL  DERIVATIVES  — 

OPRIMI  11  a  ( QIPl-OI  l/H 
OPRIMI  2  1  a  0.000 
OPRIMI S  1  •  I QIPl-OI  l/H 
APRIMI1I  a  I  AIP1-AI  l/H 
APRIMI 2  I  a  0.000 
APRIMISI  a  l  AIP1-AI  l/H 
RETURN 
END 

SUDROUTINE  C0N04I I .SHOCK .CNT ACT .J.LNOOE . RHODE .OOSTEP .RRSTEP . 
C  SSTEPI 


SUDROUTINE  CONDITION  4 

Htaaamim 

A  SITUATION  EXISTS  AT  NOOf  I  THAT  MILL  SE  CORRECTED  IN 
SUDROUTINE  CORRCT .  aNOOE  ARRAYS  ARE  6IVEN  INFORMATION 
ON  NOOE  '.OCATION  AND  SHOCK/CONTACT  SURFACE  PICTURE - 

. . VARIADLE  DEFINITIONS . 

CD4TRK  -  DENOTES  HON  MANY  NODES  NEED  TO  Of  CORRECTED 


DIMENSION  LNOOE  ( 4  )  >  RNOOE  ( 4 ) 

INTEGER  LNOOE  ,RNOOC .1 > SHOCK .CNTACT  > J .CD4TRK 
DOUBLE  PRECISION  QQSTEP ,RRSTEP  >SSTEP 

-ASSIGN  FIRST  NODE  ENCOUNTERED  DURING  THE  NEH  TIME  STEP  TO  - 

-LNOOE - 

IF( LN00EI4).LT . J )  THEN 
CD4TRX  ■  1 
END  IF 

IF( C04TRK.E9.il  THEN 
LNOOE! II  *  I 
LNOOE! 21  ■  SHOCK 
LNOOE I J  I  «  C NT ACT 
LNOOE 1 4 )  a  j 

-IF  A  SECOND  NODE  MITH  CONDITION  4  IS  ENCOUNTERED  IN  THE  - 

SHEEP,  THIS  NOOE  IS  ASSIGNED  TO  RNOOE  - 

ELSE 

RNOOE 1 1 1  •  I 
RNOOE! 2  I  a  SHOCK 
RNOOE! 5  I  a  CNTACT 
RNOOE! 4  I  a  J 

ENO  IF 

-LNOOE  ANO  RNOOE  MILL  BE  "JUMPED"  OVER  DURING  SHEEP  THUS - 

THEIR  *h*STEP  VALUES  ARE  SET  TO  0 - 

GQSTEP  ■  0.000 
RRSTEP  •  0.000 
SSTEP  ■  0.000 

IF! C04TRK.EQ.il  THEN 

CD4TRK  a  CD4TRK  ♦  I 

ELSE 

C04TRK  a  1 

ENO  IF 
RETURN 
ENO 

SUBROUTINE  C0N05I QI ,QXM1 .QIP1 ,AI .AIM! , AIP1 ,RR ,QQ,S .OELQQL , 

C  OE  LQQH , OE  LRRL . OE  LRRH . OE  LSL , OE  LSH . OE LQL . OE  LQH , 

C  OELAL  .OELAH  ,H  .EE  ,OELT  .(MINT  , PRINT  .SINT ,  AINT . 

C  QPRIH.APRIM, SHOCK, CNTACT I 


i 


SUBROUTINE  CONDITION  5 


i 


* 


i 


FOR  SUPERSONIC  FLOH  MITH  A  DISCONTINUITY  ON  ONE  SIDE  OF  THE  NOOE- 
CALCULATES  (MINT .PRINT , SINT .APRIM.QPRIM, AINT  - 

DIMENSION  LMOl  5  I.OELXI  5  I  .QINTI  3  I ,  AINT l  3  I  ,E I  3  I  .QPRIHI  3  I  .APRIMI  3  I 
INTEGER  K, SHOCK  .CNTACT 

DOUBLE  PRECISION  QI ,AI .RR ,QQ.S ,QIM1 ,QIP1 , AIM1 .AIP1 , 

C  EE  >OELT , LMO ,DELX ,9INT .AINT .E .DELQL .OELQH, OELAL 

C  QQINT . SINT, PRINT .GPRIM.APRIM, OELAH, H, 

C  DELQQL.OELQQH.OE LRRL, OE LRRH, OELSL.OE LSH 


DISCONTINUITY  ON  LEFT,  HEAOEO  RIGHT,  NOT  CROSSING  NOOE 
IF! I  SHOCK  EG. J22  I. OR.  I  CNTACT. EG. J22 l »  THEN 


INITIAL  ESTIMATE  OF  CHARACTERISTIC  SLOPES 


lmoiii  a  ex  ♦  ax 
LMDI 2 )  a  QI 
LMOIS)  ■  QIP1  -  AXP1 


vs 

A 


c 

c 

c 


- -CALCULATE  LINEARLY  INTERPOLATED  VALUES  OF  G  AND  A- 


10  K  ■  1 

20  IF  IK.LT.4)  THEN 

DELXI  K  I  a  OELT  a  LMOl K I 

eiNTIKI  *  OX  —  ( DELXI K I  »  DELQH  /  H) 

AINT(K)  a  AX  -  ( DELXI K )  »  OELAH  /  H) 

K  a  K  *  1 

SO  TO  20 

END  IF 


C 

C 

C 

c 


- CALCULATE  ERROR  BETWEEN  ESTIMATED  SLOPE  ANO  NEW  SLOPE- 

- FROM  NEW  INTERPOLATED  VALUES - 


Ell  I  «  DABS! LHOI 1 )  -  IQXNTI1I  ♦  AINTI 1 II) 

El  2)  •  OABSI LMDI 2  )  -  QINTI 2  1 1 

El  3)  ■  OABSI LMOIS)  -  I QINTI 3)  -  AINTI 3  III 


LMDI 1 1  a  QINTI 1 1  ♦  AINTI 1 1 

LM0I2)  «  QINTI 2  I 

LM0I3)  a  QINTI  3  I  -  AINTI  3  I 


C 

C 

c 


- -COMPARE  ERROR  TO  ERROR  TOLERANCE  LEVEL,  ITERATE  TO  MEET - 

IFIIEI1).6T.EE).0R.IEI2).CT.EE).0R.IEI3).ST.EEI)  GO  TO  10 


-CALCULATE  LINEARLY  INTERPOLATED  VALUES  OF  REIMAN- 
-VAR1ABLES  ANO  MODIFIED  ENTROPY  AT  POINT  A - 


QQXNT  «  QQ  -  I  DELHI  1  I  *  OELQQH  /  HI 
SINT  «S-  I0ELXI2I  *  OELSH  /  HI 
RRINT  ■  RR  -  I DELXI 3 1  •  DELRRH  /  HI 


-  CALCULATE  SPATIAL  DERIVATIVES  — 


QPRIMI1)  a  I QIP1-QI  l/H 
QPRIMI 2  I  *  0.000 
QPRIMI3I  *  I  QIPl-QI  l/H 
APRIMI 1  I  a  I AIPI-AI J/H 
APRIMI 2  I  a  0.000 
APRIMI  3  I  a  I  AIPI-AI  )/H 
END  IF 


C 

c 

c 


—  DISCONTINUITY  ON  RIGHT,  HEADED  LEFT,  NOT  CROSSING  NOOE  — 
IFI I  SHOCK. EQ. 232  I  .OR . I CNTACT . EQ . 232  I  I  THEN 
- INITIAL  ESTIMATE  OF  CHARACTERISTIC  SLOPES - 


LMDI 1  I  a  0IM1  ♦  AIM1 
LMOl  21  a  91 
LMOIS  I  •  GI  —  AI 


-CALCULATE  LINEARLY  INTERPOLATED  VALUES  OF  Q  ANO  A- 


30  K  a  1 

40  IF  IK. IT. 4  i  THEN 

DELXI K  I  a  OELT  •  LMOl K  I 

QINTI  K  I  ■  01  -  I  DELXI  K  I  a  DELQL  /  HI 

AINTIKI  a  AI  -  I0ELXIKI  a  OELAL  /  HI 

K  »  K  ♦  X 

SO  ”0  »0 

END  IF 


C 

C 


- CALCULATE  ERROR  BETWEEN  ESTIMATED  SLOPE  ANO  NEW  SLOPS  - 
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uu  u  u  u  u  ouuu  ooo  u  uuuuuuuu  o  u  uuuuuuuu 


- FROM  MEN  INTERPOLATED  VALUES - 

Ed)  ■  DABSILMOIl)  -  (QINTIll  ♦  AINTI1))) 

CIO  »  DABSILM0I2)  -  SINT  I  2  ) ) 

IIS1  a  DABSILM0I3)  -  (QINTI3)  -  AINTI 31)1 

LMDI1)  ■  QINTU)  ♦  AINTI 11 

LMDI2)  ■  QINTI  2 1 

LMD(  3  )  a  SINT  I  3 )  -  AINT131 

- COMPARE  ERROR  TO  ERROR  TOLERANCE  LEVEL,  ITERATE  TO  MEET - 

IF  1 1  El  D.GT.EE  ).  OR.  1  El  2  >.GT.  EE  >.  OR .  I  El  3  ).GT .  EE  ) )  GO  TO  30 

- CALCULATE  LINEARLY  INTERPOLATED  VALUES  OF  REIMAN - 

- VARIABLES  ANO  MODIFIED  ENTROPY  AT  POINT  A - 

MINT  •  QQ  -  I0ELXU1  *  (OELQQL  /  HI) 

SINT  «S-  I OELXI 2  I  »  I OELSL  /  Hll 
RRINT  ■  RR  -  (OELXI 3)  a  I DELRRL  /HI) 


-  CALCULATE  SPATIAL  DERIVATIVES  - 


QPRIMI 1 1 
QPRIMI 2 ) 
QPRIMI 3 1 
APRIMI 1 ) 
APRIMI 2 1 
APRIMI 3 ) 
ENO  IF 
RETURN 
END 


(QI-QIM11/H 

0.000 

IQX-QIMl )/H 
I AX-AXM1 )/H 
0.000 

I AI-AIMI  )/H 


SUBROUTINE  COND6I SHOCK ,CNTACT .HALT  I 


a 

SUBROUTINE  CONDITION  A 


INTEGER  SHOCK ,CNT ACT, HALT 

PRINT  »  , 'THE  SITUATION  FOR  CONDITION  6  KITH  THE  CONTACT  SURFACE 
PRINT  *  , 'TO  THE  RIGHT  OF  A  SHOCK  BOTH  HEADED  RIGHT  OR  THE 
PRINT  a  , 'CONTACT  SURFACE  TO  THE  LEFT  OF  A  SHOCK  BOTH  HEAOED  LEFT 
PRINT  a  , 'ADDITIONAL  LOGIC  IS  REQUIRED  TO  PROCEED* 

PRINT  a  .‘SHOCK  SHOCK,'  CONTACT  SURFACE  ■ ' .CNTACT 

MALT  a  1 

RETURN 

ENO 

SUBROUTINE  COND7I SHOCK , CNTACT ,DELT .SIGMA ,12 , 1 ,H .HALT ,Q ) 


•  • 

•  SUBROUTINE  CONDITION  7  a 

a  a 


DIWNSION  SIGMAI  4,2  1 ,121  4 ) 

INTEGER  HALT .SHOCK, CNTACT, T2, I 
DOUBLE  PRECISION  OELT .SIGMA  ,H ,9 

PRINT  a  .'CONDITION  7  REQUIRES  THAT  THE  CONTACT  SURFACE  AND 
PRINT  a  , ' SHOCK  I  MOVING  IN  OPPOSITE  DIRECTIONS  I  MEET  ANO  CROSS. 
PRINT  a  ,’NMEN  THEY  INTERSECT .THE  RESULT  IS  A  FUNCTION  OF  EXIST - 


173 


PRINT  *  ,'ING  CONDITIONS  AROUND  THEM.  DOTH  THE  ORIGINAL  SHOCK 
PRINT  »  > ' AND  CONTACT  SURFACE  HILL  EXPERIENCE  VELOCITY  CHANGES  . 
PRINT  »  , 'THIS  SUBROUTINE  MOULD  HAVE  TO  CALCULATE  HHEN  AMO  HHERE 
PRINT  ft  , 'WITHIN  THE  TIME/SPACE  ITERVAL  THE  SHOCK  AND  CONTACT 
PRINT  ft  •‘SURFACE  INTERSECT.  THIS  IS  BASED  ON  KNOW  SPEED  AND 
PRINT  ft  • 'THEIR  RESPECTIVE  LOCATIONS.  THIS  NEH  TIME •  DELT(NEH) 
PRINT  ft  •'THEN  COULD  BE  USEO  TO  RERUN  "SHEEP'S  HITH  CONDITION  7S 
PRINT  «  •'EXISTING  AT  TIME  *  T  ♦  DELTINEHI.  ADDITIONAL  LOGIC  IS 
PRINT  ft  •‘REQUIRED  TO  CONTINUE.' 

PRINT  ft  , 'SHOCK  SHOCK,'  CONTACT  SURFACE  »',CNTACT 

PRINT  ft  • 'SIGMA! 1,1 )  * • ,SIGMA< 1 ,1 ) , '  SIGMA! 1,2)  * ' , SIGMA! 1 ,2  ) 

PRINT  »  ,  'SIGMA! 2,1)  * ' .SIGMA! 2 ,1 1 , '  SIGMAI2.2)  « ' , SIGMA! 2 ,2  ) 

PRINT  »  , 'CONTACT  SURFACE  VELOCITY  *',Q 

PRINT  »  ,'NOOE  I  ft',1 

HALT  «  I 

RETURN 

ENO 

SUBROUTINE  C0N07N! SHOCK, CNT ACT ,DELT , SIGMA, 12 , 1 ,H,HALT ,Q,X ) 


ft  SUBROUTINE  CONDITION  7N  • 


ft 


DIMENSION  SIGMA! 4,2  1,12(4) 

INTEGER  HALT, SHOCK, CNT ACT, 12, I 
DOUBLE  PRECISION  DELT, SIGMA, H.Q.X 

PRINT  •  .'CONDITION  7N  REQUIRES  THAT  HHEN  EITHER  A  SHOCK  OR 
PRINT  •  .'CONTACT  SURFACE  JUMPS  A  NODE  IAS  DETERMINED  BY 
PRINT  •  .COMPARING  SIGMAIL.lt  TO  SIGMA I L  ,2  )  )  A  CONTACT  SURFACE' 
PRINT  ft  .'OR  SHOCK  MOULD  BE  MET  AND  CROSSED  DURING  THE  JUMP. 
PRINT  *  ,'TME  RESULT  HHEN  THEY  INTERSECT  IS  A  FUNCTION  OF 
PRINT  ft  .'EXISTING  CONDITIONS  AROUND  THEM.  BOTH  THE  ORIGINAL' 
PRINT  •  .'SHOCK  ANO  THE  CONTACT  SURFACE  HILL  EXPERIENCE  VELOCITY 
PRINT  •  ,  CHANGES.  THIS  SUBROUTINE  MOULD  HAVE  TO  CALCULATE  HHEN 
PRINT  ft  ,  WITHIN  THE  TIME  INTERVAL  ANO  HHERE  SPACI ALLY  THE 
PRINT  •  .'INTERSECTION  OCCURS.  THIS  IS  BASED  ON  KNOW  SPEED  ANO 
PRINT  ft  .'THE  RESPECTIVE  LOCATIONS  OF  THE  SHOCK  ANO  CONTACT 
PRINT  ft  .'SURFACE.  THE  NEH  TIME.  DELTINEHI,  COULD  THEN  BE  USED 
PRINT  •  .TO  RERUN  "SWEEP"  WITH  CONOITION  4  AT  -HIS  NODE,  ANO 
PRINT  •  >  ’ SIGMAI  2  •  2  I  *  SIGMAI1.2I  SO  THAT  CONOITION  7S  MOULD 
PRINT  •  .RESULT  IN  'HE  NEXT  TIME  INTERVAL.’ 

PRINT  ft  .SHOCK  «' .SHOCK,  CONTACT  SURFACE  «',CNTACT 

PRINT  ft  .-SIGMAIl.il  , SIGMAI  1,1  I.'  SIGMAI1.2I  ■  '  , SIGMAI  1,2) 

PRINT  •  .'SIGMAI2.lt  •' .SIGMAI  2,1 1. *  SIGMAI2.2)  •'  .SIGMAI 2 ,2  ) 

PRINT  ft  ,  'CONTACT  SURFACE  VELOCITY  *',0 

PRINT  •  ,  NOOC  I  LOCATION  AT  .<  •  '  ,X 

HALT  «  l 

RETURN 

ENO 

SUBROUTINE  CONOTS! SIGMA .HALT .SHOCK , CNT ACT  I 


SUBROUTINE  CONOITION  7S  • 

i 


DIMEWION  SIGMAI  4,2  I 
IWTE6ER  SHOCK  .CNT ACT  .MALT 
double  jre::s;dn  sigma 


PRINT  ft  , 'CONOITION  75  HAS  BEEN  MET  THIS  MEANS  THAT  THE  SMOCK 
PRINT  ■  .  ANO  CONTACT  SURFACE  ARE  LOCATEO  AT  THE  SAME  X  AT  A 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


MINT  »  ,'TIHf  OTHER  THAN  ZERO.  THIS  SUBROUTINE  MOULD  HAVE  TO  ' 
PRINT  *  , 'DETERMINE  THE  RESULT  OP  THE  INTERSECTION  RASED  ON  THE  ' 
PRINT  *  , 'CONDITIONS  TO  THE  LEFT  AND  RIGHT.  ADDITIONAL  LOGIC  IS' 
PRINT  *  , 'REQUIRED  TO  PROCEED. ' 

PRINT  •  , 'SHOCK  SHOCK,1  CONTACT  SURFACE  ■' .CNT ACT 

PRINT  m  ,'SIGMAU.l)  * '  .SIGMA  ( 1,1),'  SIGHAI2.il  •  •  ,  SIGMA  1 2 ,1 1 

HALT  •  1 

RETURN 

END 

SUBROUTINE  CONDOI SHOCK , CNT ACT , HALT  1 


SUBROUTINE  CONDITION  a  • 


INTEGER  HALT, SHOCK, CNT ACT 

PRINT  *  , ‘CONDITION  8  COULD  RESULT  ONLY  AFTER  THE  ORIGINAL  SHOCK 
PRINT  ■  ■  *  HAS  CROSSED  THE  CONTACT  SURFACE.  ANO  THE  SUBSEQUENT 
PRINT  •  , ‘CONDITIONS  DETERMINED .  ADDITIONAL  LOGIC  IS  REQUIRED  'O' 
PRINT  ■  .  CONTINUE . ' 

PRINT  •  ,‘SHOCK  ■‘.SHOCK.'  CONTACT  SURFACE  ■‘.CNT ACT 

HALT  ■  1 

RETURN 

END 

SUBROUTINE  CORRCTi  LNOOC  .RHODE  .N. SIGMA  ,H  .00  ,RR  .S  .£ .«!  .it  .12  .C  .*. 

C  AR .00 .VS. A ,9  l 


•  DISCONTINUITY  CORRfC’ION  SUBROUTINE 


.  VARIABLE  DEFINITIONS  . 

ENTRPY  -  HOC  IKIED  ENTROPY 
SNOSPO  -  SONIC  VELOCITY 

JA  -  VEL0C:,’>  RELATIVE  *0  *Hf  SHOCK,  _E»T  3JDE 
JO  -  vfLOC:’*  RELATIVE  *0  'HE  SHOCK.  RICW  SIDE 
/LC'Y  -  VELOCITY 

INTEGER  LFOOC  .RHODE  .N.  It.  NOOE  .SMOCK  ,CNT  ACT  .K 

OIHBNBION  LNOOC  I  A  I  .RNOOEi  *  I  .  SIGMA  l  A.  I  I  ,QQi  N  I  .RRt  N  I  .*1  N  I  .  1 1  •  *  >. 

C  Y?‘  •  !  .Al  N  •  .0'  N  '  .XA  •  >  .  T»i  «. 

DOUBLE  »re::sion  SICMA.JO.RR.S.a.j. 

S  H.S.Sl.GI.C  .H.AR  DQ.  .-s. 

C  RRA  .RRB  .00*  .DOB  .AA  .  AB  .QA  .OB  .SA  .  30  . 

C  SA  I  ,  SA2  .VLCTY  .SMOSRO  .ENTPPv  ,  DOC  Al  C  .  RRC  Al  C  . 

C  KO.KA 

. DC 'I  HE  STATEMENT  nPCTIONB . 

OOCALC  Vic* '  .  SNOSRO  .fNT»p.  »  nc*'  •  :NCSPt»»ENY**» 

RRC  AlC  I  Vic  *  *  .S»«SPO  ,EHT*Pv  I  •  Vic*'  SNCSKO^EMyr*. 

. SET  INI  *  I  AL  VALUES  0»  VAR1A0LES  3C  Zf  RC 

RRA  •  3  300 

RRO  •  0  000 

DBA  •  1  300 


AA  *  i  „'ij 

AO  •  0  BOO 
DA  •  0  300 


!  “  < 


*  ■  O.DOO 

SA  ■  0.000 

sa  ■  o.ooo 

00  15  Ml, 4 

XA I K  1*0.00 
XBIK  )*0.D0 

IS  CONTINUC 

- DETERMINE  IF  ONLY  ONE  NOOE  NEEDS  TO  BE  CORRECTED  OR  THO- 

- NODES.  ALSO,  IF  SHOCK  ANO  CONTACT  SURFACES  ARE  CLOSE  — 

- ENOUGH! WITH  IN  2H )  TO  INTERACT - 


IF!  I  LNOOE  I  2  l.EQ.  100  i.OR.l  LNOOE  I  5  I .  EQ .  100  I  I  THEN 

IF! I LNOOE I  1  l.GE . I RNODEI 1  1-2  I  (.AND. ( RNOOEI 1 1.GT.O  )  )  THEN 
GO  TO  20 

ELSE 

GO  TO  30 

END  IF 

ELSE  IFI RNOOE! ll.GT.OI  THEN 
GO  TO  20 

ELSE 

GO  TO  10 

END  IF 


--BRANCH  HERE  IF  ONLY  ONE  NOOE  TO  BE  CORRECTED  (LNOOE)  ANO,— 
--SMOCK/CONTACT  SURFACE  INTERACTION - 

10  SHOCK  •  LNOOE!  2  1 

CMT  ACT  *  LNOOE  <  3  ) 

-SHOCK  ON  LEFT.  HEAOEO  RIGHT,  JUMPS  OR  NOT  tCONTACT  SURFACE  ON- 
-41 GMT ,  HEAOEO  LEFT.  DOESN'T  CROSS  NODE  - 


IFI ( SMOCK. EG. 322. OR. SMOCK. EQ. 321 1. ANO.  I  CNT ACT. EG. 232  I )  THEN 


---DETERMINE  NOOE  TO  RIGHT  OF  SHOCK  AND  CONTACT  SURFACE- 

121  1  l  *  LNOOE!  11*1 
12'  2  I  >  LNOOE!  11*1 
CALL  DEL TAX  I  12 .SIGMA  .>8  ,XA  >H  I 


- EXTRAPOLATE  ”0  RIGHT  FACE  OF  CONTACT  SURFACE- 


IF  •  121  2  i.EG.Ml  then 

CALL  3B0RY I RRl N  I  ,GQ( N  )  ,S! N  )  ,RRB,QQB .SB ,AB,QB  ) 

ELSE 

CALL  EXTRAPI  RRl  12!  2  )>,RR(  12!  2  1  +  1 1.QQt  121  2  ) )  ,QQf  12!  2  1  +  1  1, 
S(  12!  2  n.S!  12!  2  )♦!  ),XB!  2  )  ,H  .RRB.QQB  ,S8  ,  AB  ,QB  ) 

ENO  IF 


- CALCULATE  VARIABLE  CHANGE  ACROSS  CONTACT  SURFACE - 

SA  «  Si  12!  2  )-I  ) 

CALL  CSJUMPt  AB . OB , SB . SA  >G2  >Q A , AA  1 

- IF  SHOCK  INTERACTS  CALCULATE  CHANGE  IN  VARIABLES  ACROSS  SHOCK - 

IF! SHOCK  EG. 321  )  THEN 

GB  *  QA 
AB  ■  AA 
SB  «  SA 

CALL  GK JUMP! AB  >38 ,38 , AR , DQ , VS ,G ,G1 .H.AA.GA ,SA I 
QQA  ■  QGCALCI OA ,AA ,SA  ) 

RRA  *  RRCALCI QA , AA ,SA  ) 


- INTERPOLATE  '0  -400EI  12!  1  l-l  I  ANO  ASSIGN  CORRECTED  7ALUES- 


IF  (  121  1  ).EQ.  2  I  THEN 

CALL  BBDRYI  RRA  ,QQA  ,SA  ,RRI  12(11-1), QQ(I2(  1  1-1  1 


St  X2(  1  )-X  )  ,At  X2t  1  )-l )  ,Q(  121 1  )-l  I ) 

ELSE 

CALL  INTERPI  RRA,RR(  X2t  1 1-2  )  ,QQA  ,QQ<  X2t  1  )-2  )  ,SA, 
St  X2<  11-2)  »XA(  ll.IXMlUHl  »RRt  12(11-1), 
QQ(  121 1  )-l ), St  X2(l)-1), At  12(11-1), 

Q(  I2(l)-1)) 

END  IF 

QQ(X2(2)-1)  ■  QQCALCI  QA,AA,SA  ) 

RRI I2( 2  )-l )  ■  RRCALCt  QA ,AA ,5A  ) 

SI  12(21-1)  ■  SA 
At  I2(  2  1-1 )  ■  AA 
3(12(21-1)  *  QA 


END  IF 


-SHOCK  ON  RIGHT,  HEADED  LEFT,  CROSSES  NODE  OR  NOT  I  CONTACT  SURFACE- 
-ON  LEFT ,HEAOED  RIGHT,  OOES'NT  CROSS  NODE - 

ELSE  IFt (SHOCK. EQ. 232. OR. SHOCK. EQ. 231). AND. (CNTACT.EQ. 322))  THEN 


-—DETERMINE  NODE  TO  RIGHT  TO  SHOCK  AND  CONTACT  SURFACE- 

1211)  *  LNOOE(l) 

12(2)  *  LNOOE(l) 

CALL  DELTA*! 12, SIGMA, X8.XA.H) 


—EXTRAPOLATE  TO  LEFT  FACE  OF  CONTACT  SURFACE- 


IF  (12(2). EQ. 2)  THEN 

CALL  BBDRYI RRI 1 ) ,QQ( 1 )  ,S( 1 )  ,RRA,QQA ,SA ,AA ,QA ) 

ELSE 

CALL  EXTRAPt  RR(  1212  )-l  5  ,RR(  I2(  2  )-2  )  ,QQ(  12(21-1)  ,QQt  12121-2), 
S(  I2(  2 1-1 1 ,S(  1 2 (  2 1-2 ) ,XA( 2  1  ,H .RRA.QQA ,SA  ,AA  ,QA  1 

END  IF 


- CALCULATE  VARIABLE  CHANGE  ACROSS  CONTACT  SURFACE - 

SA  *  SI  12(21) 

CALL  CSJUMPI  AA,QA,SA,SB,G2,QB,AB) 

—IF  SHOCK  INTERACTS  CALCULATE  CHANGE  IN  VARIABLES  ACROSS  SHOCK - 

IF( SHOCK. EQ. 231 1  THEN 
QA  a  QB 
AA  a  AB 
SA  a  SB 

CALL  GKJUMPI AA,QA,SA > AR ,DQ ,VS >G,G1 ,H,AB ,QB,S8 ) 
QQB  a  QQCALCI QB,AB, SB) 

RR8  a  RRCALCI QB , AB ,SB ) 


---INTERPOLATE  TO  NOOE(I2(D)  AND  ASSIGN  CORRECTED  VALUES- 


IF  (12(1). EQ.N)  THEN 

CALL  BB0RY(  RRB , QQB, SB, RRI  I2(  1 1 )  ,QQ(  I2(  l ) )  ,S(  12(1)) 
A(  12(1))  ,Q(  12(1))) 

ELSE 

CALL  INTERPI  RRB, RRI  I2<  1  )*1 ), QQB  ,QQ(  12(1  HI), SB, 

S<  It'  1  HI  l.XBI  1  > ,  I  XB(  1  HH  ),RR(  121  1 )  ), 

QQ(  I2(  1 1  )  ,S(  12(1))  ,A(  12(1))  ,Q(  12(D)) 

END  IF 

QQ(  I2(  2  )  I  ■  QQCALCI  08, AB, SB) 

RRI 12(21)  •  RRCALCI 06, AB, SB) 

S(  12(21  I  ■  SB 
Al  I2(  2  I  I  *  AB 
Qi  121  2  I  I  <  06 


(NO  IF 


V*OC»  ON  RIGHT. RfAOfO  RIGHT.  DOESN'T  CROSS  OR  ON  LEFT, HEADED— 


no  non  non  non  non  non 


C - RIANT,  OOCS  CROSS)  AMD  CONTACT  ON  LIFT  ,  HEADED  RIGHT ,  CROSSES--- 

C - OR  NOT!  OR  SHOCK  ON  RIGHT  >  HEADED  LIFT.  DOESN'T  CROSS  NOOC  — 

C - KITH  CONTACT  ON  LEFT,  HEAOSO  RIANT .  AND  CROSSED  NOOC  — 

C 

ELSE  If I  1 1  SHOCK . CA. ttt -OR. SHOCK . EG. 321  I . AND . <  CNT ACT . EG. 321  OR . 
C  CNT  ACT .f A. 322  )  I  .OR. I  SHOCK . EA. 212 . ANO .CNTACT . EQ. Ill  1 1 

C  THEN 

DETERMINE  NOOE  TO  RIGHT  Of  SHOCK  AND  CONTACT  SURFACE . 

1211)  •  LNOOCll)  *  1 
121  2  I  a  INODE  (  l '  ♦  l 
CALL  OELTAXI  I2.SIGHA.XB.XA.H1 

EXTRAPOLATE  TO  RIGHT  FACE  OF  SHOCK . 


IF  II2I1I.E0.N)  THEN 

CALL  BSORYI  RRIN  )  .OQI  N  >  .SI  N  >  .RRB.QGB  .SR.AG.QB  I 

ELSE 

CALL  EXT RAP!  RRI 121  1  ))  ,RRI  121  1  1*1  1 .001 121  1  )  I  .OQI  121  1  1*1  I . 

SI  121  1 1). SI  121  1 )*1  I.XBI  1  I.H.RAS.QQS.SS.AR.QR  l 

END  IF 


-CALCULATE  CHANGE  IN  VARIABLES  ACROSS  SHOCK - 


IFI SHOCK. EG. 232  I  THEN 
AA  a  AS/AR 

5A1  •  IG1/G  >*OLOG<  <  2.DOO*C*IH**2 1-0*1  000  >  /  I  6*1.  DOOM 
SA2  ■  G1*OLOGI  I  (G-l. 000  l»l  H**2  1*2.000  >  /  110*1. 000  >• 

I H**2 I  I  I 
SA  a  S8«SA1*SA2 
US  •  06  -  VS 
UA  a  ue  -AAaOQ 
QA  a  UA  ♦  VS 


END  IF 


CALL  SXJUMPIAB.QB.SS.AR.DO.VS.G.Gl.H.AA.QA.SA) 


-CALCULATE  CHANGE  IN  VARIABLES  ACROSS  CONTACT  SURFACE  - 


IFI CNTACT .EO. 322 )  THEN 

OQI  12(1 1-1 )  a  OQCALCIOA.AA.SA  > 
RRI 121 1>-1)  a  RRCALCtOA.AA.3A) 
SI  121  11-1  I  *  SA 
Al  121  1  1-1  I  a  AA 
QII2ID-1)  a  QA 


08  a  QA 
AB  a  AA 
SB  a  SA 

SA  a  S)  121  2  1-2  I 

CALL  CSJUMPI AB.QB.S8.SA.G2.QA.AA) 
QQA  a  QQCALCIQA.AA.SA) 

RRA  a  RRCALCf QA.AA.5A  ) 


-INTERPOLATE  TO  NOOEI 121 2  )-l  1  ANO  ASSIGN  CORRECTED  VALUES- 


IF  (12(2). EQ. 2)  THEN 

CALL  BBORYI  RRA  .QQA  ,SA  ,RRI  121  2  )-l )  ,001 121  21-1), 
SI  121  2  )-l ) ,  Al  121  2  )-l )  ,Q(  121  21-1)1 

ELSE 

CALL  INTERPI  RRA , RRI  121  2  )-2  I.QQA.OOI  121  2  )-2  I.SA, 

SI  121  2  )-2  )  iXAl  2  )  ,1 XAI  2  )*H  )  ,RR(  121  2  1-1 ), 

QQI 121  2  )-l )  ,SI  121  2  )-l )  ,Al  121  2  )-l ) , 

0(  I2(  2  )-l ) ) 

ENO  IF 


END  IF 


—SHOCK  ON  LEFT,  HEADEO  LEFT,  DOESN'T  CROSS  NOOE,  OR  ON  RIGHT,— 


ono  non  non  o  o  n 


-m wo  l«*t,  cnuii  MOM >  ttm  contact  on  rxomt,  mcaoco  left,— 

-CKtUI  MOM  OM  NOT:  00  SMOCK  ON  LEFT,  MS  A0€0  RI8MT .  DOESN'T  — 
-CMOS  MOOt .  ANO  CONTACT  ON  RJONT .  Hi  AMD  LETT.  CHOOSE  0  NOOC  — 


ELSE  XFI  M SMOCK  . C«. lit  . OK. SMOCK. CO. ESI  I . AMO . (CNTACT . 10. 212  OR . 
C  CNTACT.C0.2S1  II  .00 . 1  SMOCK . 10. 122 . AMO .CNTACT . (0. 2J1  I  I 
C  TNCN 


— MTCRHXM  MOM  TO  RXOMT  Of  SMOCK  ANO  CONTACT  SURFACE  - 

X2I 1  I  «  LNOOC I  X  I 
121  2  I  •  LNOOC  l  1  l 
CALL  OtLTAXi  I2.SIGMA.>«.XA.MI 


-EXTRAPOLATE  ’0  LEFT  FACE  OF  SMOCK  - 


XF  l  121  1  I. CO. 2  I  TMCN 

CALL  BOORV  I  Ml  1  1 .001  1  I  .SI  1  I  ,MA  .00 A  .SA  . AA  ,0A  I 

ELSE 

CALL  EXTRAS!  RRI  121  1  1-1  I.AAl  X  21  1  1-2  1.001  121  1  1-1  1 .001  121  1  1-2  I. 
SI  121  1  1-1  I.Sl  121  1  1-2 I.XAI  1  I.H.MA.QOA.SA.AA.QA  I 

CNO  XF 


--CALCULATE  CHANCE  IN  VARIABLES  ACROSS  SMOCK  - 


IFI  SMOCK.  CO.  122  I  TMCN 
AO  ■  AA/AR 

SA1  a  (  01/0  JR0L06I  I  2  .  DOO**C*l  Mm2  1-0*1  ■  000  I  /  (6*1.000  I 
SA2  ■  61*0L06I  I  10-1.000  I»IM*mi2  1*2.300  i  /  IIG*'.D00I* 

I  M*m*Z  I  I  • 

SO  a  SA»SA1«SA2 

UA  a  01  •  VS 

UO  a  IIA  -  18*00 

CO  a  UO  *  VS 

ELSE 

CALL  SKJUNRI  AA.OA.SA.AR.OO.VS.S.Ol.M.AO.OO.SOI 

CNO  IF 


—CALCULATE  CHANCE  IN  VARIABLES  ACROSS  CONTACT  SURFACE - 

XFl  CNTACT  .  EQ.  232  I  THEN 

QQIX2I  11  I  •  OQCALCIOO.AB.S8  I 
RRI  X2I  1  I  I  ■  RRCALC10B.AB.SBI 
SI  12'  l  l  I  a  SB 
Al  121  1  I  I  *  IB 
01  X2I  1  M  >  ge 


AA  <  AS 
SA  a  SB 

SB  <  Si  121  2  1*1  ) 

CALL  CSJUMP1AA,OA,3A,SB,G2,GO,A8I 
006  a  QQCALCIQB.AB.SB  I 
RRB  >  RRCALCIQB.A8.S6  I 

-INTERPOLATE  TO  NOOCI 12(21)  ANO  ASSIGN  CORRECTED  VALUES- 


EHO  XF 


IF  I  121  2  I.CO.NI  THEN 

CALL  880RYIRRB.QQ6.S6.RRI  121  2  I  I  >001 121  2  I  i  .Si  I2i  2  I  I , 
Al  121  2  )  I  >Q<  1212)11 

ELSE 

CALL  INTERPt  RRB. RRI  121  2  1*1  l  ,006.001  Id  2  1*1  1,58, 

SI  121  2  1*1  I.XBI  2  I.IXBI  2  i*H  I, RRI  121  2  I  I , 

00(121  2  1 1, SI  121  2  I),  Al  121  2  I  1,01 121  2  1 1  I 

END  IF 


m 


£ 


; 


» 


r*  ■■ 
VJ 


C - MUNCH  HERE  IF  LNOOE  AND  RNOOE  AM  CLOU  ENOUGH  FOR  SHOCK - 

C  - AND  CONTACT  SURFACE  INTERACTION - 

C 

C  LIFT  NOOCi 

C  ---SHOCK  ON  RIGHT,  HEADED  RIGHT,  JUMPS  NOOE  KITH  CONTACT  ON  LEFT— - 

C  - HEADED  RIGHT,  JUMPS  NOOE  OR  NOT |  OR  NO  SHOCK  HITH  CONTACT  ON  - 

C  —LEFT,  HEAOED  RIGHT,  JUMPS  NOOE  — 

C  RIGHT  NOOE: 

C  —SHOCK  ON  LEFT,  HEAOEO  RIGHT,  JUMPS  NOOE  HITH  NO  CONTACT  SURFACE  — 

C 

20  IF)  ( I  LNOOEl  2  I  .EQ.  221 ) .  ANO.  (  LNOOEI  I). EG. 121. OR.  LNOOEI I ).  EG.  122  )  I 
C  .OR. t  LNOOEI  2  I  EQ. 100 . ANO . LNOOEI  S I . EQ. S21 1 »  THEN 

XFI I RNOOEI 2  I.EQ. J21  I . AND . I RNOOE  I  5 I.EQ. 100 l I  THEN 
C 

C - DETERMINE  NOOE  TO  RIGHT  OF  SHOCK  ANO  CONTACT  SURFACE . 

C 

121 II  ■  RNOOEI II  *  1 
IFILNOOEI2l.EQ.S22i  THEN 
121 2  I  ■  LNOOEI 1 1 

ELSE 

121  2  I  ■  LNOOEI 1 1  ♦  1 

ENO  IF 

CALL  QELTAXI X2  ,SIOIA  ,XB  .XA  ,H  ) 

C 

C  - EXTRAPOLATE  TO  RIGHT  FACE  OF  SHOCK . 

C 

IF  <12111. EQ. HI  THEN 

CALL  MORY I  RRI  N  I  ,QQ(  N  I  ,S<  N  )  ,RR&  ,QQB  >S8 ,  AB  ,06  ) 

ELSE 

CALL  EXTRAPI  RRI  121  1  I  I  .RRI  121  l  )♦!  1 .001  121  1  i  1 .001  121  1  1*1  I, 

C  S<  121  1  I  I, SI  121  1  1*1 I.ABI  1  >, H.RRB, 906. SB, AB. 06  I 

ENO  IF 
C 

C  . CALCULATE  VARIABLE  CHANCE  ACROSS  SHOCK . 

C 

CALL  SKjUMPI AS,Q6.SB.AR,0q,VS,G,Gl,H,AA,9A,SA) 

C 

C . DETERMINE  CORRECTED  VALUES  AT  NOOEl 121 1  )-l I  ANO  ASSIGN  THEM - 

C 

IFI LNOOEI 2  I.EQ. 100  I  THEN 

QQA  •  OQCALCt  QA . AA.3A  I 
RRA  *  RRCALCIQA.AA.SA  I 

C 

C . INTERPOLATE  to  NOOE I  121 1  I- 1  I  ANO  ASSIGN  CORRECTED  VALUES - 

C 

IF  (  121  1  I.EQ. 2  I  THEN 

CALL  6B0RYI  RRA  ,QQA  ,SA  ,RRI  121  1  1-1  ),QQI  121  1  )-l  I, 

C  SI  121  1 1-1  I, At  121  1  )-l  )  ,01  121  1  )-l  )  I 

ELSE 

CALL  INTERPI  RRA, RRI  121  1  1-2  l.OOA.QQl  121  1  1-2  I,  SA  , 

C  SI  121  1  i-2  I.XAI  l  I.IXAI  1  i*M  l, RRI  121  1  l-l  I, 

C  QQl  121  1 1-1  I, SI  121  1  (-1 1.AI  121 1  l-l  I, 

C  Q(  12(1  l-l  H 

ENO  IF 
C 

C  . EXTRAPOLATE  TO  RIGHT  FACE  OF  CONTACT  SURFACE - 

C 

IF  I  121  2  I.EQ.NI  THEN 

CALL  BBORYi  RRI  N I ,QQ(N I ,S(N  I ,RRB ,006 ,S8,AB,Q6  I 

ELSE 

CALL  EXTRAPI  RRI  121  2  I  )  ,RRI  121  2  1*1  I, QQl  I2l2n,QQII2l2l  +  ll, 
C  31  121  2  I  1,31  121  2  !♦!  I  ,XBi  2  I  ,H  ,RRB  ,9QB,S6  ,AB  ,9B  l 


ENO  IF 

QQl  121  11-11  *  9QCALCI  QA  ,AA  ,SA  ) 
RRI  12(11-11  *  RRCALCIQA.AA.SA  I 
SI  121  11-11  *  SA 
Al  121  1  1-1)  ■  AA 
Ql  121 11-1)  ■  QA 


180 


nnn  ono  non  nnnnnn  nnn 


r 


XPiLNooiisi.io.m)  tmcn 

90(12(2  )>  ■  «0(I2l  11-1 1 
M(  12(2)1  ■  Ml  X2(  1  )-l ) 

SI  12(2)1  ■  SI  121 1  )-l ) 

*112(2)1  ■  *112(11-1) 

0112(2))  ■  01 121  D-D 
00  TO  40 

I  LSI 

0*  ■  0* 

*0  ■  A* 

SO  ■  SA 

l NO  IP 

EN0  IP 

S*  »  S(  12(21-2  I 

CALL  CSJUHPI AS. 00. SO. SA >02.0*. A* I 
00A  ■  QOCALCIQA.AA.SA) 

HR*  •  RRCALCI QA,AA ,SA  ) 

INTERPOLATE  TO  N00E( 121 2 1-1 )  AND  ASSIGN  CORRECTEO  VALUES - 

IP  (X2(  2  1.EQ.2)  THEN 

CALL  BOORYI  RRA  ,0QA  ,SA  ,RR(  1212  1-1 1 ,001 12(  21-11, 
C  SI  121  2  )-l  >,AI  121  2  1-1 1, 0(  121  21-1)1 

ELSE 

CALL  1NTERPI  RRA  ,RR|  121  2  1-2  )  .00*  ,QQI  1212  1-2  I  .SA , 

C  SI  121  2  1-2  >,XAI  2  I.IXAI  2  )*H  I  ,RR(  121  21-11, 

C  001 12(2 1-1 1,31 12(2  1-11, At  121  2  1-11, 

C  01 121 2  )-l )  I 

END  IP 

ENO  IP 


LEFT  NOOE i 

- SHOCK  ON  LEFT,  HEAOEO  RIGHT,  JUMPS  NOOE  HITH  NO  CONTACT- 

RIGHT  NOOE: 

—NO  SHOCK.  HITH  CONTACT  ON  RIGHT,  HEAOEO  LEFT, JUMPS  NOOE  — 

ELSE  IF( ( LNOOEI 2  ).EQ. 321 1 . ANO . ( LNOOEI S  1  .EO. 100  1 1  THEN 

IPI ( RNOOEI 2  1 . 29. 100  I  .AND. I RNOOEI J  1  .EQ. 231 1 1  THEN 

- DETERMINE  NOOE  TO  RIGHT  OP  SMOCK  ANO  CONTACT  SURFACE . 

12(1)  «  LNOOEI  1)  ♦  1 
1212)  «  RNOOEI 1) 

CALL  OELTAXI 12, SIGMA, H8, HA, H  I 

- CALCULATE  JUMP  THROUGH  CONTACT  SURFACE  THEN  SHOCK - 

IPI  121 1  ).E0<(  121  2  1-1 1 1  THEN 
AA  »  A<  121  1  I  I 
QA  ■  Ql  121  1  I  I 
SA  a  SI  121  1  I  I 
SO  ■  SI  121  2  1*1 1 

CALL  CSJUMPI AA ,QA,SA ,S8 ,G2 ,00 ,AS  1 
000  •  QQCALCI 00 i AS , SO  1 
RRS  ■  RRCALCI  00, AS, SO  I 

- INTERPOLATE  TO  NOOEt 12(211  ANO  ASSIGN  CORRECTEO  VALUES - 

IP  II2l2l.Eq.NI  THEN 

CALL  BOORYI  RRS  .006  ,SB  ,RRI  12(2)1 ,901 121  2  1 1  ,SI  121  2  1 1 , 
C  Al  12(21)  ,QI  12(21)1 

ELSE 

CALL  INTERPI  RRO.RRI  121  2  1*1  1 ,900,001 121  2  1*1  1  ,SB, 

C  SI  121  2  )♦!  ),XBI  2  1,1X01  2  ItHl.RRI  121  2  1 1, 

C  OQI 121  2  1 1, SI  121  2  1  l.AI  121  2  11,01 121  2  1 )  1 

EM)  IF 

AO  ■  Al  121  1 1 ) 

08  a  0(1211)1 
SB  ■  SI  121  1)1 
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.AR, 09. VS. 9. 91  .M.AA.9A.SA  I 


99A  •  99CALd9A.AA.SA  I 
MU  ■  R9CAICI9A.AA.SAI 


-INTIRPOLATI  TO  MOO«l  IK  1 1-1 1  AND  AS3I0N  CORRECTEO  VAUCS- 


ir  iimi.ii.t)  thkn 

CAU  990RV<RRA.99A,SA,RRII(I  1 1-1 1. Ml  IK  11-11. 
«l  IK  11-11  .Al  IK  11-11 .01  IK  1 1-1 II 

•  LSI 

CALL  INTIRPIRRA.RRU2I  1 1-2  I.MA.MI  IK  1  I-2I.SA. 
SI  IK  1  1-2 I.XAt  l  I.IXAI  1  l»M  I.RRI  121  1  »-l  I. 

<M<  121 1 1-1  I. Si  121  1  1-1  i. At  I2i  1  l-l  i, 

91 121  1  l-l  1  I 

(NO  IF 


Ml  121 11-1 1  •  RRI 121 1 1-2  I 
Ml  IK  11-1 1  •  Ml  IK  11-2  l 
SI  IK  11-1 1  ■  SI  121 11-2  l 
Al  121 11-1  I  •  Al  121  11-2  l 
•I  121  11-1 1  •  91 121 11-2  I 
Ml  121  2  1 1  •  RRI  121  2  1*1  I 
Ml  121  21  I  •  Ml  IK  2  1*1  I 
SI  121  2  1 1  ■  Si  121  2  >*1  I 
Al  121  2  II  •  Al  I2l  2  1*1  I 
91 121  2  1 1  •  911212  1*11 


(NO  If 


(NO  IF 


RIOHT  NOOI : 

—  SHOCK  ON  RIOHT,  HEAOEO  LIFT.  JUMPS  NOOI  KITH  NO  CONTACT  — 
LIFT  NOOI : 

—NO  SMOCK.  NITM  CONTACT  ON  LIFT.  HCAOIO  RIOHT.  JUMPS  NOOI-  — 


(LSI  IFMRNOOK  2  I.I9.2SI  I.ANO.IRNOOK  1  1.(9. 100  i  I  THIN 

IF I  I LNOOII 2  I  ■ (9. 100  I . ANO . I L NOOI l 1 1 . 19. S21  I  I  THIN 


-OETIRMINE  NOOI  TO  RIOHT  TO  SMOCK  ANO  CONTACT  SURFACE  - 


121  1  I  •  ANOOE I  1  I 

12(21  •  LNOOII II  ♦  1 

CALL  OILTAXI 12. SIGMA. IM.XA.H I 


-CALCULATE  JUMP  THROUGH  CONTACT  SURFACE  THIN  SHOCK  - 


IFII2I2  i.EO.i  IK  1  l-l  i  I  THIN 
AS  ■  Al  121  2  1 1 
09  ■  91  IK  III 
S9  •  SI  IK  2  i  I 
SA  •  SI  121  2  t-2  I 

CALL  CSJUMPi  A9.Q0.39.SA.02.9A.AA  I 
99A  ■  9QCALCI9A.AA.SA  I 
RRA  ■  RRCALCI9A.AA.SA  I 


-INTERPOLATE  TO  NOOK  IK  2  l-l I  ANO  ASSIGN  CORRECTED  VALUIS- 


IF  112121.(9.2)  THIN 

CALL  SOORYI  HR  A,  GO  A,  S  A,  RRI  121  2  l-l  1 ,901  1 21  2  l-l  I , 
SI  IK  2  l-l  I  .Al  121  2  l-l  i.Qi  IK  2  1-1)1 

ELSE 

CALL  IMTERPI  RRA, RRI  121  2  1-2  I.OOA.OQI  121  2  S-2  I.SA, 

SI  IK  2  1-2  I  ./Al  2  )  ,i  XAl  2  l«H  i  .RRl  I2<  2  l-l  l . 

Ml  121  2  l-l  I, SI  121  2  l-l  I , Al  121  2  l-l  I, 

9112121-1)1 

(NO  IF 

AA  *  Al  121  2  I  I 
OA  •  91  121  2  I  I 
SA  ■  SI  121  2  I  I 

CALL  SKJUMPI  AA.9A,SA,AR,09.VS.O,C1,M,A9,M.S9I 


o  u  u  g  u  u  u  u  u  u  u  u  u  g  u 


too  •  oocalcioo.ao.ooi 

«•  ■  99CAUI00.A9.90I 

C 

c - nnwouti  to  nook  hi  hi  am  usiw  coomctio  valuoo 

c 


If  I  III  II. It. Ml  TMtN 

CALL  0009VI990. 900. 90. Ml  HI  lll.99<12<  111. SI  121 II  I. 
C  Al  III  II  1,01  Xtl  I  Ml 

ILM 

CALL  INTI 991  MO  .Ml  1(1  1 1*11.000.001  III  1 1*1 1. SO. 

C  01X11  1 1*1  I. KOI  1 1. 1  Ml  1  UN  I, Ml  III  1  I  I. 

C  001  ltl  1  > I.Sl  121  1 1  I.AII2I  1 11.01  121  1  l  l  I 

(NO  It 

(LSI 

Ml  12t  2  l-l  i  •  Ml  121  2  1-2  l 
OOl  221  2i-li  ■  Mi  I2<  2  1-2  t 


SI  ltl  2  1-1  l  ■  SI  121  2  1-2  l 
Al  121  2  1-11  •  Al  121  2  1-2  i 
91 121  2  l-l  l  •  91  121  2  1-2  i 
Ml  121  1  1 1  *  Ml  12!  1  l«l  l 
OOl  Z 21  1  l  I  •  OOl  III  1  HI  I 
Si  III  1  1 1  *  Si  111  1  ul  l 

Al  121  1 II  •  Al  111  1  l«l  i 

9l  121  1  I  l  ■  91 121  lull 

(NO  If 

(NO  If 
C 

C  010MT  NOOfl i 

C  ---SMOCK  ON  HfT .  MIND  L«fT .  JUNTA  NON  M1TN  CONTACT  ON  RIQMT--- 

C  ---•KAMO  LIFT,  dots  NOT  CIIOSS  NOM  39  NO  SMOCK  MUM  CONTACT  ON--- 

C  ---A1QMT .  rtf  AMO  LIFT.  CKOS  HO  fON  — 

C  LIFT  mom 

C  ---SMOCK  ON  KI9MT.  rtfAOlO  LIFT.  Jim  NOM  MUM  NO  CONTACT--- 

C 


ILSf  IFi  I  iNNOOIl  2  '.(9. Ill  '  ANO.lRNOM' 1  1(9.211  09  RNOOU  1  I  (9. 
C  212H  09  iNNOMi  2  1.(9. 100  ANO  KNOM  •  1  l  .  (Q .  211  i  l  TIKN 

If  i  i  LNOMi  2  I  19  211  I  .ANO  ILNOOU  D  IO.  100  i  l  THIN 


MTI9MN0  NOM  TO  9I0MT  Of  SMOCK  AND  CONTACT  SUN  FAC  I 

111  1  I  •  LNOMI  1  I 
I F  l  KNOM  i  1  i  19  2I«l  TNCN 

III  2  i  •  KNOMi  ll*l 

(LSI 

III  2  l  •  KNOMi  1  i 

(NO  IF 

CALL  MLTAXI  It.OIONA.jK.XA.rt  I 
(XTKAKOLATI  TO  LIFT  FACI  Of  SMOCK - 


If  i  III  V  i  (9. 2  i  tm«N 

CALL  0009YIN9I  1  l.9Q<  1  l. SI  1  1 .99  A  .90A.SA  .AA.9A  I 

ILM 

CALL  IXT9A9I  Ml  X 21  1  l-l  1.991  121  1  1-2  1 ,90t  IHli-1  1.901  III  1  l-l  I, 
C  SI  121  1  l-l  I.SI  121  1  1-2  I.XAI  1  I  ,M  ,RKA  >9QA  >SA  .  AA  ,  9A  I 

(NO  If 


CALCULATI  /ARIAOLI  C HANOI  ACROSS  SMOCK . 

CALL  SKJUN9I AA,9A.SA.A9,09./S.S.01.M,A0.90,SO  I 

■  Of  T( RHINO  COKRICTfO  VALUIS  AT  N0MU2I1H  ANO  ASSIGN  TM(M 

IFI9N0MI2I.I9.1001  THIN 

009  •  0QCALCI99.A0.S9  I 
MO  •  ARC  ALC I  09  ■  AO  •  SO  i 

IMTI990LAT!  TO  N00III2I1II  ANO  ASSX6N  C099ICTIO  VALUIS - 
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n  n  r>  nonoo  r»oo  non 


It  (Itl  11. IS. Ml  TMN 

CALL  RSORYl  RM  ,0M  ,M  .Ml  Itl  X  11  .Ml  Itl  1 1  I.SI  Xtl  111. 
Aixtuii.sixtnin 

|  CM 

CALL  IMTIMI  Ml. Ml  Itl  1  Ull.Mt.MUtl  1 1«1 1.M. 

SI  Xtl  1 1*1 1  .Ml  1  I  >1  Ml  lUMI.MIXtilll. 

SSI  Xtl  1 1). SI  Xtl  1 1  l,AI  Xtl  1 1 ). Si  Xtl  1 II I 

■mo  it 


c 

c 

c 


-IXTRAROLATI  TO  LIPT  FACt  OF  CONTACT  SURFACE- 


It  IItltl.tS.il  TNCN 

CALL  MOST  I  Ml  1  I  .Ml  1  I.SI  1  I  .AAA  .MA  .SA  .AA  .SA  I 

I  LSI 

CALL  (XTRARiRRl  Itl  1 1-1  I  .RSI  Itl  t  i-t l.QQi  Itl  t  1-1 1. SSI  Itl  t  l-t  I. 
SI  Itl  t  1-1  I.SI  Xtl  t  l-i  I.XAI  t  I.H.MA.SQA.SA.AA.4A  l 

(NO  It 
(LSI 

Ml  Xtl  1 1 1  •  SSCALCIM.AS.SS  I 
Ml  Itl  II  I  •  MCALCl  SS.AS.SS  I 
SI  Itl  HI  •  SS 
Al  Itl  1  I  I  •  AS 
SI  Itl  1  I  I  •  OS 
Itl RNOOf  < S  I . (S . til  I  thin 

•SI  Xtl  1  1*1 1  ■  Ml  Itl  1  II 
Ml  Itl  1  1*1  1  ■  SQI  Xtl  1  I  I 

si iti n«i i  •  sum ii 

Al  Xtl  1  I*  1  I  •  Alltllll 
SI  Itl  11*1  1  •  Ollttlll 
SO  TO  <*0 

■  LSI 


SA  •  SS 


I  NO  It 


■  NO  It 


SS  •  Si  Xtl  1 1*1 1 

— CALCULATI  VARXASLI  CMANSI  ACROSS  CONTACT  SURFACE- 

CALL  CSJUHRI  AA.SA.SA.SS.Ot.M.AI 
SM  •  OQCALCiaS.AS.SSI 
•M  •  MCALCl  OB. AS. SS* 


-  INTI  MOL  ATI  '0  NOOK  xtl  tn  AND  ASSIGN  CORMCTIO  VALUC3- 


It  I  Itl  t  I.IS.MI  THIN 

CALL  SSORVIRRS.SM.SS.RRI  Itl  1 1 1. Ml  Xtl  1 1 1  .SI  Itl  1 1 1 . 
Antitii.ontitiii 
(LSI 

CALL  IMTIMI  MS. Ml  Itl  t  i*l  l  .90S >901  Itl  t  1*1  I  .SS. 

SI  Itl  t  1*1  I.XSl  t  I.IXBI  t  l*M  I  .Ml  Itl  t  II. 

Ml  Itl  1 1  I  .SI  IK  t  I  I  .Al  Itl  t  I  1. 01  Itl  1 1  1 1 

(NO  It 


(NO  It 

■no  xr 

00  TO  40 

—  BRANCH  HIRf  It  THt SI  ARC  ONI  OR  THO  NOOIS  TO  CORRICT .  BUT  THIY  — 

—  ARC  SIM  RATIO  ST  MORI  THAN  XH.  CHICK  FOR  SHOCK/CONTACT  SURFACI  — 

—  I  NT  I A  ACT  ION,  ANO  CORRICT  AMRORRI ATILY  — 

SO  NOOI  ■  LNOOI I  1  I 

SMOCK  ■  LNOOIItl 
C NT ACT  «  LNOOI 1 1 1 

—  NO  SMOCK,  CONTACT  SURFACI  ON  LIFT,  HI AO I 0  RIGHT,  JUMPS  NOOI  — 

SS  Xtll SMOCK .10- 100  I  .AND. I CNTACT .10. Stl  1 1  THEN 
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c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


zt(ti  •  moo e  ♦  x 

CALL  0CLTAX<X2.SX«HA.>«,XA.H> 

—cmcx  FOR  A  SHOCK  INTERACTING  MXTH  CONTACT  SURFACE  AT  NOOC - 

XFM  SXOHAIX  .S > .  LT .  I  Xlt  2  )*H  II .  AND . 

C  I$2GNA(1.2).9T.(X2I2>>)>TH(N 

AS  a  Al 1212)1 

or  ■  91x21211 
SR  a  SI  12(21 1 

(LSI 

— EXTRAPOLATE  TO  RIGHT  FACE  OF  CONTACT  SURFACE - 

XF  ( 121 2  1 .10. Nt  THIN 

CALL  SRORY t  RRI N ) ,091  N  I  ,S<  N I  .RRR  .990  .SR .  AR  .OR  I 

(LSI 

CALL  (XTRAPI  RRI  X2<  2  1 I  .RRI  X2I  2  I«1 1 .90*  12(2  11  .991  121  2  (♦!  I  • 
C  SI  121 2  )  I.SI  III  2  1*1  I.XRI  2  I.H.RRR.GGR.SR.AR.QR  I 

(NO  IF 

(NO  IF 

SA  a  Si  12(2  1-2  1 

-  CALCULATI  VARXARLl  CHAN6I  OVIR  CONTACT  SURFACI  — - 

CALL  CSJUHPI  AR.9R.SR.SA  ,62  .9A.AA ) 

9QA  a  RQCALCI QA  ,AA  >SA ) 

RRA  a  RRCALCIQA.AA.SA  ) 

--INTERPOLATE  TO  NOOK  121  2 l-l )  ANO  ASSIGN  CORRICTEO  VALUES - 

IF  1 121  21. IQ.  2  I  THIN 

CALL  SRORYI RRA.QQA .SA .RRI 121 2  l-I  I  .991 121 2  1-1  I . 

C  SI  121  2  1-1 1.AI  121  2  »-l  I.QI  X2I  2  >-l  >  > 

ELSE 

CALL  INTIRPI  RRA  .RRI  X2I  2  »-2  I  >QQA  ,QQI  121  2  1-2  )  .SA , 

C  $<  121  2  >-2  I  .XAI  2<. IXAI2I*HI.RR  1X21  21-1). 

C  QQI  III  2  l-l  I  .SI  121  2  1-1 1  .At  III  21-11, 

C  Q(  X2I  2  l-l  1 1 

(NO  IF 

-  NO  SHOCK.  CONTACT  surface  on  lift,  headed  LEFT.  JUMPS  nooe  — 

(LSI  IF l l SHOCK . E9. 100  ) . ANO ■ I CNTACT .19.211)1  THEN 
121  2  I  ■  NOOE 

CALL  OELTAXI  1 2  .SIGMA  ,)•  .XA  ,H  I 

- CHECK  FOR  A  SHOCK  INTERACTING  HITH  CONTACT  SURFACE  AT  NOOE . 

IF)  l  SIGMA  l  1.2  I  ST.  1X21  2  l-l  2  .  DOO«H  I  I  I.ANO.  I  SIGMA)  1.2  I.LT. 

C  1X21 21-HIII  THEN 

AA  ■  Al  121  2  1—1 ) 

QA  ■  91 121  21-11 
SA  •  SI  121  2  l-l  I 

(LSI 

- EXTRAPOLATE  TO  LEFT  FACE  OF  CONTACT  SURFACE . 

IF  112(21.(9.2)  THEN 

CALL  ROOKY (  RRI X  1 ,9Q< 1 )  ,Sl 1 )  .RRA .QQA.SA.AA.QA  ) 

ELSE 

CALL  EXTRAPIRRI  121  21-1 1, RRI  121  2  1-21,991 12*2  »-l  » .00*  12'  2  >-2  » . 

C  SI  X2I  2 )  —  1 >,S(  121  2 1  —  2  I ,XA I  2  > ,H ,RRA .QQA.SA.AA.QA  I 

END  IF 
I  NO  IF 

SB  »  Si  121  2  1*1  i 

CALL  CSJUMPt  AA,QA,SA,SB,G2,QR,AB) 

QQQ  a  QQCALCI  QB.AB.SB  ) 

RRR  a  RRCALCI  QB.AB.SB) 
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I 


cs 


S 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


- INTERPOLATE  TO  NOOE  I  IK  2  ) )  ANO  ASSIGN  CORRECTED  VALUES - 

IP  IIKII. 10. Ml  THIN 

CALL  MONVI  KM  >000  .SO. Ml  Hill  1 .00*  IK  1 1)  .SI  IK  l  II , 

C  AII2t2)I.G<XK2lll 

■  LSI 

CALL  INTERPt  RM.RRI  III  1 1*1 1 .000.001  IK  X  Ml  i.SO. 

C  S(I2ltl*ll.XDI2>.l»l2MMI.RRII2l2H. 

C  001  III  X  1 1  .SI  121  X  II  .Al  121  2  II  .01 121  2  1 1 1 

END  IF 

-  SMOCK  ON  LEFT.  MEAOCO  KXCMT .  JUMPS  NOOE  KITH  NO  CONTACT  — 

ELSE  XFItSMOCK.EO.S2ll.ANO.ICNTACT.EO.100n  THEN 

X  2 1  1 1  •  NOOE  «  1 

CALL  DELTAXI  12, SIGMA  .M.XA.H  I 

- -CHECK  FOK  A  CONTACT  SUKFACS  INTERACT INO  MITM  SHOCK  AT  NOOE - 

XFI  I  SXCMAI  2 .2  I .  LT  .  I  X21 1  MM  I  I .  ANO . 

C  ISIGMAI  2,2  I.GT  .1X21  1  I  I  I  I  THEN 

AO  •  Al  121  1 1  1 
00  ■  01X21  II  I 
SO  ■  SI  X  21 1 1 1 

ELSE 

- EXTRAPOLATE  TO  RIGHT  FACE  OF  SMOCK - 

IF  tXZIll.EO.Nl  THEN 

CALL  OOORYI RRI N 1 ,001 N I  .SI  N I  ,RRO  .000  .SO  .AO  .00  I 

ELSE 

CALL  EXTRAPI  RRI  121 1 1  I.RRI  121 1 1*1  l.OQI  X21 1  >  1 .001 12 1  1  Ml  I , 
C  SII2Illl,Sll2tl>*ll,XBll).H.RRB.GQO,SO.AS.QOI 

ENO  IF 

ENO  IF 

—  CALCULATE  VARIAOLE  CHANGE  ACROSS  SHOCK  --- 

CALL  SKJUMPI  AS.OO.SO.AR.OO.VS.G.Gl.H.AA.OA.SAI 
QQA  •  OOCALCIOA.AA.SA  I 
HR  A  *  RRCALCIQA.AA.SA  I 

- INTERPOLATE  TO  NOOE I  121 1 *- 1  I  ANO  ASSIGN  CORRECTED  VALUES . 

IF  II2lll.EO.2l  THEN 

CALL  OOORYI RRA.OOA.SA, RRI  121 1 1-1 1 ,001 121 1 1-1 1 , 

C  SI  X21 1 1-1 1  .Al  121 1 1-1 1  ,Qt  121 1 1-1 1 1 

ELSE 

CALL  INTERPI  RRA.RRI  121  1  1-2  I  .QQA.OQI  I2(  1  1-2  I  ,SA, 

C  StI2lll-2I.XAIll,IXAIlMHl,RRII2<ll-ll, 

C  OQI 121 1 1-1 1  ,SI  121 1 1-1  I.AI  121 11-11, 

C  01 121 11-1 1 1 

ENO  IF 

—  SHOCK  ON  RIGHT,  HEADED  LEFT.  JUMPS  NOOE  KITH  NO  CONTACT  SURFACE— 

ELSE  XFI I  SHOCK . EQ. 2S1 1 . ANO . I CNTACT . EQ . 100  I  I  THEN 
121  1 1  ■  NOOE 

CALL  OELTAXI 12 .SIGMA ,XB  >XA ,H  I 

IF  I  ISIGMAI  2,2  I.GT.IX21 1  1-1 2 . DOO»H  1 1 1 . ANO . I SIGMAI  2,2  I.  LT. 

C  I  X2I  1  l-H  I  I  I  THEN 

AA  ■  Al  121  1 1-1  I 
QA  ■  0II2IH-11 
SA  »  SI  121  1  1-1  I 

ELSE 

- EXTRAPOLATE  TO  LEFT  FACE  OF  CONTACT  SURFACE - 


C 

c 

c 
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UUU  u  u  u  uuu  uuu  u  uuuuuuu  uuuuuu  uuo 


e 


SA.AA.QA) 


xf  iieid.eg.si  tncn 

CALL  INRYI  Ml  1 1  ,QQ(  1 )  ,SI  1 )  ,RRA  >Mt  , 


ELSE 

CALL  CXTRAPI  MUM  )-l  1  iMI  X21 1  l-l  1 ,Q0< 121 1 1-1 1 ,991  X21 1 1-2  >. 


SI  X21 1 1-1 1 ,51  X21 1 1-2  I  ,XA(  1 1 ,H,RIU ,QQA,3A ,AA ,QA  ) 

END  IP 
(NO  XP 


-  CALCULATt  VARIABLE  CHANGE  OVER  SMOCK  - 


XPI  X21 1 )  .IQ.N-1 1  THEN 
I  NO  XP 

CALL  SKJUMPI  AA.QA.SA.AR.OG.VS.a.Gl.H.AB.QB.SBI 


-OCTBRHINC  CORRECTED  VALUES  AT  N00EIX2ID)  ANO  ASSIGN  THEN' 


IPIX2I1I.EQ.N-1)  TNCN 
END  XP 

QOS  ■  QQCALCI OS, AS, SO ) 

RRS  ■  RRCALCI  QB.AS.SS  ) 

- INTERPOLATE  TO  NOOEII2I1II  AND  ASSIGN  CORRECTED  VALUES - 

IP  1 121 1 ) .EQ.N I  THEN 

CALL  SSORYI  RRS  ,999  ,S8  ,RRI  X21 1 )  )  ,QQ(  X21 1 ))  ,SI  X21 1 1 1 , 
C  At  X21 1 1 )  ,QI  X2I  lilt 

ELSE 

CALL  XNTERPI  RRS  ,RRt  X2I  D«l>  .QQB.QQI  121 1  )*1 )  ,SS , 

C  SIX21 1  Ull.XBI  11.IXBI  1  )«H ),RRI  121 1) ), 

C  QQI  12(111  ,S)  121 1 1 )  ,A(  X21 1  I )  .Ql  12(D)) 

CNO  XP 

END  IF 

-  CHECK  IF  SECOND  NODE  NEEDS  CORRECTING  IP  SO  LOOP  BACK  AROUND  — 

XPIRNOOCID.EQ.O)  GO  TO  AO 
NODE  «  RNOOE(l) 

SHOCK  •  RNOOEI  2  ) 

CNTACT  •  RNOOEI D 
RNOOEI 11  •  0 
GO  TO  IS 

AO  RETURN 
ENO 

SUBROUTINE  DCLTAXI 12 .SIGMA ,XB,XA  ,H  ) 


*  LOCATE  DISCONTINUITY  HITMIN  INTERVAL  SUBROUTINE  * 


DIMENSION  121 A ) , SIGMA I A ,2  ) ,XSI A ) ,XAI 4  ) ,X2I A ) 

INTEGER  12 

OOUBLE  PRECISION  SIGMA, >®,XA,H,X2 
X21 3  1*0.000 
X2I A  1*0.000 

- CALCULATE  DISTANCE  NOOE  TO  RIGHT  OP  DISCONTINUITY  IS  FROM  LEFT.  — 

- BOUNDARY  OF  TUBE  THEN  DETERMINE  DISTANCE  FROM  THIS  NOOE  121*1 - 

- TO  DISCONTINUITY  ANO  OISTANCE  FROM  NOOE  TO  LEFT  OF  DISCONTINUITY- 

- TO  THAT  DISCONTINUITY - 

X2I1)  ■  0SLE<X2(D-D  *  H 
XBU)  »  i  X2<  1 1-SIGMAI  1,2)1 
XAI  1  I  *  (H-XBI  D  I 

- ^CALCULATE  SAME  VALUES  FOR  CONTACT  SURFACE - 


MUUUUUUUUUUUUUU  u  u u  uou  u  u  u  u  uouoouu  ouuu  o  ouuu 


Xtttl  •  0BLEH2<21-I>  •  N 

xbixi  •  ix«2t-siam<2,<n 
xtiti  •  iH-aitn 
RETURN 
(NO 

SUBROUTINE  SKJUMPt  AON, 00M, SOM. ARATIO, DCLTAR, VSHOCK , 6.01  .N, 
C  AUP,  GUP.  SOOt 


•  SMOCK  JUMP  SUBROUTINE 


»•***•***■ 


- VARIABLE  DEFINITIONS - 

■ON  -  VARIABLE  OOFGttTREAM  OF  SHOCK 
ARATIO  -  SPEED  OF  SOUND  RATIO 
OELTAG  -  VELOCITY  JUMP  ACROSS  SMOCK 
■UP  -  VARIABLE  UPSTREAM  OF  SHOCK 
VSHOCK  -  SHOCK  VELOCITY 

DOUBLE  PRECISION  AON .QON .SON, ARATIO, DELTA* .VSHOCK ,6 .61 ,H, 
C  AUP ,QUP »SUP (SA1 ,SA2 .UON  >UUP 

- -CALCULATE  THE  SPEEO  OF  SOUND  CHANCE  THRU  SHOCK - 

AUP  ■  AONoARATIO 

- CALCULATE  THE  ENTROPY  CHANGE  THRU  SHOCK - 

SA1  ■  (  Gl/Q  l>OLOGI  (  2 . 0O0«G»l  H**2  l-G+l .  DOO  )/(  G*1 .  DOO  )  I 
SAX  ■  G1*OLOG(  (  ( 6-1 . 000  N**2  >♦2  1/11  6*1 . 000  M*»2  )  1  ) 

SUP  »  SON  -  SA1  -  SA2 

- -CALCULATE  THE  VELOCITY  CHANCE  ACROSS  SHOCK - 

UON  «  QON  -  VSHOCK 

UUP  ■  UON  ♦  A0N*0ELTA« 

CUP  *  UUP  ♦  VSHOCK 

RETURN 

ENO 

SUBROUTINE  CSJUMP( AON, QON, SON, SUP ,S2 ,QUP ,AUP > 


*  CONTACT  SURFACE  JUMP  SUBROUTINE  * 


■  ■ 


DOUBLE  PRECISION  QUP  » QON. AUP , AON > SUP  ,SDN,G2 

- CALCULATE  THE  SPEEO  OF  SOUND  CHANGE  ACROSS  THE  CONTACT  SURFACE  — 

- NOTE!  THE  VELOCITY  IS  CONTINUOUS  ACROSS  A  C.S. . 

CUP  ■  QON 

AUP  «  AON  ■  DEXPf I  SON-SUP  1/621 

RETURN 

ENO 

SUBROUTINE  EXTRAPI  RRL .RRR  >QQL ,QQR ,SL ,SR ,DX ,DH , PRINT >QQINT ,SINT , 
C  AINT.QINT) 


■ 

■ 

■ 


EXTRAPOLATION  SUBROUTINE 


■ 

■ 

■ 


non  oooonor*  n  o  r>  n  o  nnoooo 


C 


C 

c 

c 

c 


c 


MUKi  PRECISION  RRL  .RRR .RRL  .RM.SL  .SR .  ax .  OH.  RRIMT  .MINT  .SINT  . 

C  41  NT  .RIKT 

- -REIMAN  VMX4M.CS  4M  EXTRAPOLATED,  AFC  TMi  CORRESPONDING . 

- VALUSS  FOR  VELOCITY  4MB  MID  OF  KUC  4 M  CALCULATED . 

MINT  •  ML  -  I  OX  •  I  MM -ML  1/OM  I 
MINT  •  ML  -  (OX  ■  I  MR -ML  I/M  ! 

SINT  ■  St  -  I DX  •  I SR-SL  I/DM  I 
4XMT  «  I  MINT  -RRIMT  1/1  2 . 000*SXMT  I 
•IMT  •  I  MINT  *RRIMT  l/l  1 . 000  I 
Rt TURN 
(NO 

SUBROUTINE  INTI  Ml  URL  .RRR  .ML  .MR  .SL  .SR  .OX.DH  .RRIMT  .MINT  ,SINT  . 
C  AIMT.GIMTl 


INTERPOLATION  SUBROUTINE 


DOUBLE  PRECISION  RRL  .RRR  .ML  .MR .SL  .SR , OX, DM, RRIMT  .MINT  .SINT , 
C  AXNT.QXNT 

REIMAN  VARIABLES  ARE  INTERPOLATED.  AND  THEIR  CORRESPONDING - 

VALUES  OF  VELOCITY  ANO  SPEEO  OF  SOUND  ARE  CALCULATED - 

RRIMT  a  RRL  -  I  OX  •  i RRL-RRR  l/OH  I 
MINT  ■  ML  -  (OX  •  <  ML -MR  l/OH  I 
SINT  a  SL  -  (DX  •  I SL-SR  l/OH  I 
AINT  a  IMINT-RRINT  I/I  2.000*SINT  1 
BINT  a  l  MINT+RRINT  1/2.000 
RETURN 
ENO 

SUBROUTINE  OBURSTl N,H,M.RR,S,G,G1 ,62 ,0ELT ,12 ,X2 ,N,AR .DQ.VS, 
RLMPRES.SIGMA.A.Q) 


a  * 


«  DIAPHRAM  BURSTING  SUBROUTINE  » 


INTEGER  N.X  >Y ,12 ,L .SHXOIR .CSDIR.LHPRES 

DIMENSION  X21 4 1  ,RR<  N I  ,M(  N )  ,SI  N 1 ,12(4 ) , SIGMAI 4, 2  )  ,AI  N I  ,QI  N ) 

DOUBLE  PRECISION  X2 ,X.H . AB  >SA .SB  > AA ,QA  .06  >MA , QOS ,RRA ,RRB .SIGMA , 
C  RR  >M  ,S ,  TS  ,H  ,00 ,  AR  ,PR  ,G  ,GI  ,G2 .  VS  ,0ELT  ,CSRMN .  A , 

C  Q .MREIMN.DREMN.EREIMN.HM.SAl ,SA2 .SAP ,S8P ,XB ,XA 

♦♦♦♦♦♦♦♦  LOCATING  THE  NODE  TO  THE  RIGHT  *********** 

DO  10  L*1 .4 

SIGMAI  L.1)*SIGMAIL. 2) 

Y*0 

X=H 

1*2 

11  IF  < .NOT.IY.EQ.OI)  GOTO  10 

IF  ISIGMAlL.ll.LT. XI  THEN 
X2(L)=X 
I2ILI-I 
Y=1 
ENO  IF 
X=X*H 
I»I»l 
GOTO  11 
10  CONTINUE 


189 


noon  oonoo  o  onnnonnn 


„♦.♦♦♦♦♦♦♦♦♦♦  CORRECT  FOR  DIAPHRAM  BURSTING  ♦♦♦♦♦♦♦♦♦♦♦♦♦ 


- AT  TINE  ZERO  DETERMINE  CORRECT  SHOCK  DIRECTION - 

- SHF DIR  »  J  IS  A  SHOCK  HEADED  LEFT,  ANO  SHKDIR  3  2  IS  SHOCK - 

- TRAVELING  RIGHT - 

IF( LHPRES.EQ.3  )  THEN 
SHKDIR  «  3 
X2(l>  3  X2(l)  -  H 
12(1)  *  12(1)  -  1 
X2(  2  )  3  X2<  2  I  -  H 
12(2)  3  12(2)  -  1 
GO  TO  20 

ELSE 

SHKDIR  *  2 

EM)  IF 

20  RRA«RR(  12(1 1-1) 

RRB*RR( 121 1 )  I 
QQAzQQI 1 2 ( 1  1-1) 

QQB*Q9(  12(D) 

SA*S( 12(11-1) 

S83S(  I2(  1 H 

21  AB*( QQB-RRB )/( 2.000*S8 ) 

AA«( QQA-RRA )/( 2 . D00*SA  ) 

IF( SHKDIR. EQ. 3  )  THEN 

MREIMN  3  ( RRB-RRA  )/AA 
OREMN  3  DABS! MREIMN) 

ELSE 

MREIW  3  <  QQA-QQB  l/AB 
DREMN  3  MREIW 

END  IF 

- ITERATE  FOR  PROPER  VALUE  OF  H  USING  THE  QUADRATIC  FIT  OF  THE - 

- REIMAN  VARIABLE  CHANGE  WITH  V  CURVE.  NOTE  LEFT  MOVING  SHOCKS— 

- ARE  USED  IN  THESE  EQUATIONS  SINCE  RRB-RRA/AA=-< QQA-QQB/ AB  ) - 

100  HH=t(  3.0396408D01-(  (  OREMN+2 . 7574000  I/O.  286337D00  ) ) 
H=5.513294D00-boQRT(WW) 

0Q=2 . 000*( W*W-1 . 000  )/( H*( G+l . 000  ) ) 

AR=0SQRT(  2 . 000*1  G-l .  DOO  1*1  1.000+1  < G-l .000  )*W*W/2 .000  1  )* 

C  (G*G2*H»H-1.000  I  )/(  (G+l. 000  l*H  ) 

PR=(  2.000*G/(  G+l. DOO  )  )*W*H-(  (G-l.  000  )/IG+l.uOO  )  ) 

0Rs(  (G-l. 000  )*W*H+2.D00  )/(  (G+1.000  )»W*H) 

EREIMN=0Q+( AR-1. 000  )*G2-( AR*G1/G  )*DLOG( PR»I  DR**G  ) ) 

IF  (OABS(EREIMN-DABS( MREIMN)).LT. 0.10-5)  GO  TO  110 
DREW  3  (OABSI  MREIMN)  -  EREIMN)  +  OREMN 
GOTO  100 

- DETERMINE  S  BEHIND  SHOCK  ANO  THEN  CALCULATE  THE  RIEMAW . 

- VARIABLE  CHANGE  ACROSS  THE  CONTACT  SURFACE - 

110  SA1  3  (G1/G)*OLOG((2.DOO*G*(H**2)-G+1.DOO)/(G+1.DOO)  i 

SA2  3  Gl*DLOGI  (  (  G-l .  000  )*(  H**2  1  +  2  )/(  >  G+l .  DOO  )•(  N»*2  I  )  I 
IF( SHKDIR. EQ. 2  I  THEN 

SAP  3  SB  -  SA1  -  SA2 
S8P  3  SAP 

ELSE 

S8P  3  SA  -  SA1  -  SA2 
SAP  3  SBP 

.  END  IF 

103  IFI SHKDIR. EQ. 2)  THEN 

CSRMN  3 1  I  OEXP(  (  S8P-SA  '/G2  '  >*'  SA  '  -  S8P 

ELSE 

CSRW  «(  I  OEXP)  (  SAP -SB  i/G2  t  i*i  SB  -  :**•  > 

END  IF 

IF( SHKDIR. EQ. 3  I  THEN 


>  I 


